{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "sys.path.append('/home/yang_liu/python_workspace/3DGS')\n",
    "\n",
    "import yaml\n",
    "import torch\n",
    "import torch_scatter\n",
    "import torchvision\n",
    "import wandb\n",
    "import time\n",
    "import inspect\n",
    "import imageio\n",
    "import numpy as np\n",
    "import open3d as o3d\n",
    "import pynvml\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "from arguments import GroupParams\n",
    "from scene import LargeScene\n",
    "from os import makedirs\n",
    "from gaussian_renderer import render_lod, render\n",
    "from utils.general_utils import safe_state\n",
    "from utils.large_utils import which_block, block_filtering\n",
    "from utils.sh_utils import SH2RGB\n",
    "from argparse import ArgumentParser\n",
    "from arguments import ModelParams, PipelineParams, get_combined_args\n",
    "from gaussian_renderer import GaussianModel\n",
    "from torch.utils.data import DataLoader\n",
    "from utils.camera_utils import loadCamV2\n",
    "from transforms3d.euler import mat2euler, euler2mat\n",
    "from transforms3d.quaternions import mat2quat, quat2mat\n",
    "from utils.general_utils import PILtoTorch\n",
    "from scene.cameras import Camera\n",
    "from scipy.spatial import ConvexHull\n",
    "from matplotlib.backends.backend_agg import FigureCanvasAgg\n",
    "\n",
    "WARNED = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BlockedGaussian:\n",
    "\n",
    "    gaussians : GaussianModel\n",
    "\n",
    "    def __init__(self, gaussians, lp, range=[0, 1], scale=1.0, compute_cov3D_python=False):\n",
    "        self.cell_corners = []\n",
    "        self.xyz = None\n",
    "        self.feats = None\n",
    "        self.max_sh_degree = lp.sh_degree\n",
    "        self.device = gaussians.get_xyz.device\n",
    "        self.compute_cov3D_python = compute_cov3D_python\n",
    "        self.cell_ids = torch.zeros(gaussians.get_opacity.shape[0], dtype=torch.long, device=self.device)\n",
    "        self.mask = torch.zeros(gaussians.get_opacity.shape[0], dtype=torch.bool, device=self.device)\n",
    "\n",
    "        self.block_dim = lp.block_dim\n",
    "        self.num_cell = lp.block_dim[0] * lp.block_dim[1] * lp.block_dim[2]\n",
    "        self.aabb = lp.aabb\n",
    "        self.scale = scale\n",
    "        self.range = range\n",
    "\n",
    "        self.cell_divider(gaussians)\n",
    "        self.cell_corners = torch.stack(self.cell_corners, dim=0)\n",
    "\n",
    "    def cell_divider(self, gaussians, n=4):\n",
    "        with torch.no_grad():\n",
    "            if self.compute_cov3D_python:\n",
    "                geometry = gaussians.get_covariance(self.scale).to(self.device)\n",
    "            else:\n",
    "                geometry = torch.cat([gaussians.get_scaling,\n",
    "                                      gaussians.get_rotation], dim=1)\n",
    "            self.xyz = gaussians.get_xyz\n",
    "            self.feats = torch.cat([gaussians.get_opacity,  \n",
    "                                    gaussians.get_features.reshape(geometry.shape[0], -1),\n",
    "                                    geometry], dim=1).half()\n",
    "            \n",
    "            for cell_idx in range(self.num_cell):\n",
    "                cell_mask = block_filtering(cell_idx, self.xyz, self.aabb, self.block_dim, self.scale)\n",
    "                self.cell_ids[cell_mask] = cell_idx\n",
    "                # MAD to eliminate influence of outsiders\n",
    "                xyz_median = torch.median(self.xyz[cell_mask], dim=0)[0]\n",
    "                delta_median = torch.median(torch.abs(self.xyz[cell_mask] - xyz_median), dim=0)[0]\n",
    "                xyz_min = xyz_median - n * delta_median\n",
    "                xyz_min = torch.max(xyz_min, torch.min(self.xyz[cell_mask], dim=0)[0])\n",
    "                xyz_max = xyz_median + n * delta_median\n",
    "                xyz_max = torch.min(xyz_max, torch.max(self.xyz[cell_mask], dim=0)[0])\n",
    "                corners = torch.tensor([[xyz_min[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_max[2]]], device=self.xyz.device)\n",
    "                self.cell_corners.append(corners)\n",
    "    \n",
    "    def get_feats(self, indices, distances):\n",
    "        out_xyz = torch.tensor([], device=self.device, dtype=self.xyz.dtype)\n",
    "        out_feats = torch.tensor([], device=self.device, dtype=self.feats.dtype)\n",
    "        block_mask = (distances >= self.range[0]) & (distances < self.range[1])\n",
    "        if block_mask.sum() > 0:\n",
    "            self.mask = torch.isin(self.cell_ids, indices[block_mask].to(self.device))\n",
    "            out_xyz = self.xyz[self.mask]\n",
    "            out_feats = self.feats[self.mask]\n",
    "        return out_xyz, out_feats\n",
    "\n",
    "    def get_feats_ptwise(self, viewpoint_cam):\n",
    "        out_xyz = torch.tensor([], device=self.device, dtype=self.xyz.dtype)\n",
    "        out_feats = torch.tensor([], device=self.device, dtype=self.feats.dtype)\n",
    "\n",
    "        homo_xyz = torch.cat([self.xyz, torch.ones_like(self.xyz[..., [0]])], dim=-1)\n",
    "        cam_center = viewpoint_cam.camera_center\n",
    "        viewmatrix = viewpoint_cam.world_view_transform\n",
    "        xyz_cam = homo_xyz @ viewmatrix\n",
    "        self.mask = (xyz_cam[..., 2] > 0.2)\n",
    "        if self.mask.sum() == 0:\n",
    "            return out_xyz, out_feats\n",
    "\n",
    "        distances = torch.norm(self.xyz - cam_center[None, :3], dim=-1)\n",
    "        self.mask &= (distances >= self.range[0]) & (distances < self.range[1])\n",
    "        if self.mask.sum() > 0:\n",
    "            out_xyz = self.xyz[self.mask]\n",
    "            out_feats = self.feats[self.mask]\n",
    "        return out_xyz, out_feats\n",
    "\n",
    "def load_gaussians(cfg, config_name, iteration=30_000, load_vq=False, device='cuda', source_path='data/matrix_city/aerial/test/block_all_test'):\n",
    "    \n",
    "    lp, op, pp = parse_cfg(cfg)\n",
    "    setattr(lp, 'config_path', cfg)\n",
    "    lp.source_path = source_path\n",
    "    lp.model_path = os.path.join(\"../output/\", config_name)\n",
    "\n",
    "    modules = __import__('scene')\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        if 'apply_voxelize' in lp.model_config['kwargs'].keys():\n",
    "            lp.model_config['kwargs']['apply_voxelize'] = False\n",
    "        gaussians = getattr(modules, lp.model_config['name'])(lp.sh_degree, device=device, **lp.model_config['kwargs'])\n",
    "        scene = LargeScene(lp, gaussians, load_iteration=iteration, load_vq=load_vq, shuffle=False)\n",
    "        print(f'Init {config_name} with {len(gaussians.get_opacity)} points\\n')\n",
    "\n",
    "    return gaussians, scene\n",
    "\n",
    "def parse_cfg(cfg):\n",
    "    lp = GroupParams()\n",
    "    op = GroupParams()\n",
    "    pp = GroupParams()\n",
    "\n",
    "    for arg in cfg['model_params'].items():\n",
    "        setattr(lp, arg[0], arg[1])\n",
    "    \n",
    "    for arg in cfg['optim_params'].items():\n",
    "        setattr(op, arg[0], arg[1]) \n",
    "\n",
    "    for arg in cfg['pipeline_params'].items():\n",
    "        setattr(pp, arg[0], arg[1])\n",
    "    \n",
    "    return lp, op, pp\n",
    "\n",
    "def loadCamV3(args, id, cam_info, resolution_scale, xyz=None, z_dim=None, yaw=None):\n",
    "    # use appointed pitch and height\n",
    "    orig_w, orig_h = cam_info.image.size\n",
    "\n",
    "    if args.resolution in [1, 2, 4, 8]:\n",
    "        resolution = round(orig_w/(resolution_scale * args.resolution)), round(orig_h/(resolution_scale * args.resolution))\n",
    "    else:  # should be a type that converts to float\n",
    "        if args.resolution == -1:\n",
    "            if orig_w > 1600:\n",
    "                global WARNED\n",
    "                if not WARNED:\n",
    "                    print(\"[ INFO ] Encountered quite large input images (>1.6K pixels width), rescaling to 1.6K.\\n \"\n",
    "                        \"If this is not desired, please explicitly specify '--resolution/-r' as 1\")\n",
    "                    WARNED = True\n",
    "                global_down = orig_w / 1600\n",
    "            else:\n",
    "                global_down = 1\n",
    "        else:\n",
    "            global_down = orig_w / args.resolution\n",
    "\n",
    "        scale = float(global_down) * float(resolution_scale)\n",
    "        resolution = (int(orig_w / scale), int(orig_h / scale))\n",
    "\n",
    "    resized_image_rgb = PILtoTorch(cam_info.image, resolution)\n",
    "\n",
    "    gt_image = resized_image_rgb[:3, ...]\n",
    "    loaded_mask = None\n",
    "\n",
    "    if resized_image_rgb.shape[1] == 4:\n",
    "        loaded_mask = resized_image_rgb[3:4, ...]\n",
    "\n",
    "    if yaw is not None and z_dim is not None and xyz is not None:\n",
    "        Rt = np.zeros((4, 4))\n",
    "        Rt[:3, :3] = cam_info.R.transpose()\n",
    "        Rt[:3, 3] = cam_info.T\n",
    "        Rt[3, 3] = 1.0\n",
    "\n",
    "        C2W = np.linalg.inv(Rt)\n",
    "        euler = np.array(mat2euler(C2W[:3, :3]))\n",
    "        euler[z_dim] = yaw\n",
    "        C2W[:3, :3] = euler2mat(*euler)\n",
    "        C2W[:3, 3] = xyz\n",
    "        Rt = np.linalg.inv(C2W)\n",
    "\n",
    "        R = Rt[:3, :3].transpose()\n",
    "        T = Rt[:3, 3]\n",
    "    else:\n",
    "        R = cam_info.R\n",
    "        T = cam_info.T\n",
    "\n",
    "    return Camera(colmap_id=cam_info.uid, R=R, T=T, \n",
    "                  FoVx=cam_info.FovX, FoVy=cam_info.FovY, \n",
    "                  image=gt_image, gt_alpha_mask=loaded_mask,\n",
    "                  image_name=cam_info.image_name, uid=id, data_device=args.data_device)\n",
    "\n",
    "def loadCamV4(args, id, cam_info, resolution_scale, xyz=None, angle=None):\n",
    "    # use appointed pitch and height\n",
    "    orig_w, orig_h = cam_info.image.size\n",
    "\n",
    "    if args.resolution in [1, 2, 4, 8]:\n",
    "        resolution = round(orig_w/(resolution_scale * args.resolution)), round(orig_h/(resolution_scale * args.resolution))\n",
    "    else:  # should be a type that converts to float\n",
    "        if args.resolution == -1:\n",
    "            if orig_w > 1600:\n",
    "                global WARNED\n",
    "                if not WARNED:\n",
    "                    print(\"[ INFO ] Encountered quite large input images (>1.6K pixels width), rescaling to 1.6K.\\n \"\n",
    "                        \"If this is not desired, please explicitly specify '--resolution/-r' as 1\")\n",
    "                    WARNED = True\n",
    "                global_down = orig_w / 1600\n",
    "            else:\n",
    "                global_down = 1\n",
    "        else:\n",
    "            global_down = orig_w / args.resolution\n",
    "\n",
    "        scale = float(global_down) * float(resolution_scale)\n",
    "        resolution = (int(orig_w / scale), int(orig_h / scale))\n",
    "\n",
    "    resized_image_rgb = PILtoTorch(cam_info.image, resolution)\n",
    "\n",
    "    gt_image = resized_image_rgb[:3, ...]\n",
    "    loaded_mask = None\n",
    "\n",
    "    if resized_image_rgb.shape[1] == 4:\n",
    "        loaded_mask = resized_image_rgb[3:4, ...]\n",
    "\n",
    "    if angle is not None and xyz is not None:\n",
    "        Rt = np.zeros((4, 4))\n",
    "        Rt[:3, :3] = cam_info.R.transpose()\n",
    "        Rt[:3, 3] = cam_info.T\n",
    "        Rt[3, 3] = 1.0\n",
    "\n",
    "        C2W = np.linalg.inv(Rt)\n",
    "        euler = np.array(mat2euler(C2W[:3, :3]))\n",
    "        euler = angle\n",
    "        C2W[:3, :3] = euler2mat(*euler)\n",
    "        C2W[:3, 3] = xyz\n",
    "        Rt = np.linalg.inv(C2W)\n",
    "\n",
    "        R = Rt[:3, :3].transpose()\n",
    "        T = Rt[:3, 3]\n",
    "    else:\n",
    "        R = cam_info.R\n",
    "        T = cam_info.T\n",
    "\n",
    "    return Camera(colmap_id=cam_info.uid, R=R, T=T, \n",
    "                  FoVx=cam_info.FovX, FoVy=cam_info.FovY, \n",
    "                  image=gt_image, gt_alpha_mask=loaded_mask,\n",
    "                  image_name=cam_info.image_name, uid=id, data_device=args.data_device)\n",
    "\n",
    "# find the a & b points\n",
    "def get_bezier_coef(points):\n",
    "    # since the formulas work given that we have n+1 points\n",
    "    # then n must be this:\n",
    "    n = len(points) - 1\n",
    "\n",
    "    # build coefficents matrix\n",
    "    C = 4 * np.identity(n)\n",
    "    np.fill_diagonal(C[1:], 1)\n",
    "    np.fill_diagonal(C[:, 1:], 1)\n",
    "    C[0, 0] = 2\n",
    "    C[n - 1, n - 1] = 7\n",
    "    C[n - 1, n - 2] = 2\n",
    "\n",
    "    # build points vector\n",
    "    P = [2 * (2 * points[i] + points[i + 1]) for i in range(n)]\n",
    "    P[0] = points[0] + 2 * points[1]\n",
    "    P[n - 1] = 8 * points[n - 1] + points[n]\n",
    "\n",
    "    # solve system, find a & b\n",
    "    A = np.linalg.solve(C, P)\n",
    "    B = [0] * n\n",
    "    for i in range(n - 1):\n",
    "        B[i] = 2 * points[i + 1] - A[i + 1]\n",
    "    B[n - 1] = (A[n - 1] + points[n]) / 2\n",
    "\n",
    "    return A, B\n",
    "\n",
    "# returns the general Bezier cubic formula given 4 control points\n",
    "def get_cubic(a, b, c, d):\n",
    "    return lambda t: np.power(1 - t, 3) * a + 3 * np.power(1 - t, 2) * t * b + 3 * (1 - t) * np.power(t, 2) * c + np.power(t, 3) * d\n",
    "\n",
    "# return one cubic curve for each consecutive points\n",
    "def get_bezier_cubic(points):\n",
    "    A, B = get_bezier_coef(points)\n",
    "    return [\n",
    "        get_cubic(points[i], A[i], B[i], points[i + 1])\n",
    "        for i in range(len(points) - 1)\n",
    "    ]\n",
    "\n",
    "# evalute each cubic curve on the range [0, 1] sliced in n points\n",
    "def evaluate_bezier(points, n):\n",
    "    curves = get_bezier_cubic(points)\n",
    "    seg_length = [np.linalg.norm(points[i] - points[i + 1]) for i in range(len(points) - 1)]\n",
    "    mean_length = np.mean(np.array(seg_length))\n",
    "    return np.array([fun(t) for length, fun in zip(seg_length, curves) for t in np.linspace(0, 1, max(int(n * length / mean_length), 1))])\n",
    "\n",
    "def evaluate_bezier_v2(points, n):\n",
    "    curves = get_bezier_cubic(points)\n",
    "    seg_length = [np.linalg.norm(points[i] - points[i + 1]) for i in range(len(points) - 1)]\n",
    "    mean_length = np.mean(np.array(seg_length))\n",
    "    seg_pts_num = np.array([max(int(n * length / mean_length), 1) for length in seg_length])\n",
    "    acc_pts_num = np.cumsum(seg_pts_num)\n",
    "    pts = np.array([fun(t) for length, fun in zip(seg_length, curves) for t in np.linspace(0, 1, max(int(n * length / mean_length), 1))])\n",
    "    return pts, acc_pts_num\n",
    "\n",
    "def evaluate_bezier_v3(points, n):\n",
    "    curves = get_bezier_cubic(points)\n",
    "    seg_samples = [n]\n",
    "    speed_start = None\n",
    "    speed_end_pre = curves[0](1) - curves[0](1-1/n)\n",
    "    for idx in range(1, len(curves)):\n",
    "        speed_start = curves[idx](1/n) - curves[idx](0)\n",
    "        n = int(n * np.linalg.norm(speed_start) / np.linalg.norm(speed_end_pre))\n",
    "        seg_samples.append(n)\n",
    "        speed_end_pre = curves[idx](1) - curves[idx](1-1/n)\n",
    "        \n",
    "    acc_pts_num = np.cumsum(np.array(seg_samples))\n",
    "    pts = np.array([fun(t) for samples, fun in zip(seg_samples, curves) for t in np.linspace(0, 1, max(samples, 1))])\n",
    "    return pts, acc_pts_num\n",
    "\n",
    "def np_move_avg(a,n,mode=\"same\"):\n",
    "    b = a.copy()\n",
    "    b[n//2:-n//2] = np.convolve(a, np.ones((n,))/n, mode=mode)[n//2:-n//2]\n",
    "    return b\n",
    "\n",
    "def np_move_avg_v2(a,n,mode=\"same\"):\n",
    "    b = a.copy()\n",
    "    if len(b.shape) == 1:\n",
    "        b[n//2:-n//2] = np.convolve(a, np.ones((n,))/n, mode=mode)[n//2:-n//2]\n",
    "    else:\n",
    "        for i in range(b.shape[1]):\n",
    "            b[n//2:-n//2, i] = np.convolve(a[:, i], np.ones((n,))/n, mode=mode)[n//2:-n//2]\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 No LoD, MatrixCity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading camera 741/741\n"
     ]
    }
   ],
   "source": [
    "load_vq = True\n",
    "iteration = 30_000\n",
    "custom_test = '../data/matrix_city/aerial/test/block_all_test'\n",
    "if load_vq:\n",
    "    iteration = None\n",
    "\n",
    "config = \"../config/block_mc_aerial_block_all_lr_c36_loss_5_50_lr64_vq.yaml\"\n",
    "model_path = os.path.join('../output', os.path.basename(config).split('.')[0])\n",
    "with open(config) as f:\n",
    "    cfg = yaml.load(f, Loader=yaml.FullLoader)\n",
    "    lp, op, pp = parse_cfg(cfg)\n",
    "    setattr(lp, 'config_path', config)\n",
    "    if lp.model_path == '':\n",
    "        lp.model_path = model_path\n",
    "\n",
    "with torch.no_grad():\n",
    "    modules = __import__('scene')\n",
    "    model_config = lp.model_config\n",
    "    gaussians = getattr(modules, model_config['name'])(lp.sh_degree, **model_config['kwargs'])\n",
    "\n",
    "    if custom_test:\n",
    "        lp.source_path = custom_test\n",
    "        filename = os.path.basename(lp.source_path)\n",
    "    scene = LargeScene(lp, gaussians, load_iteration=iteration, load_vq=load_vq, shuffle=False)\n",
    "\n",
    "    bg_color = [1,1,1] if lp.white_background else [0, 0, 0]\n",
    "    background = torch.tensor(bg_color, dtype=torch.float32, device=\"cuda\")\n",
    "    \n",
    "    if custom_test:\n",
    "        views = scene.getTrainCameras() + scene.getTestCameras()\n",
    "\n",
    "org_scaling = gaussians._scaling.clone()\n",
    "org_opacity = gaussians._opacity.clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cam_center_list = []\n",
    "dim = 1\n",
    "for idx in range(len(views)):\n",
    "    Rt = np.zeros((4, 4))\n",
    "    Rt[:3, :3] = views[idx].R.transpose()\n",
    "    Rt[:3, 3] = views[idx].T\n",
    "    Rt[3, 3] = 1.0\n",
    "\n",
    "    C2W = np.linalg.inv(Rt)\n",
    "    cam_center_list.append(C2W[:3, -1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.1 Get Convex Hull of Training / Test Poses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f94f4df8b50>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAH0CAYAAACzc26fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7p0lEQVR4nO2deXjU1PrH365QRBEoW6ELBcq+C4giFFBARUUR4SpcFbef1+t+Ra6KylURFTdwQUC5V1xQARdkuSJcRJaisoogBYFOBWpZZeneeX9/pJlkpkkmmUkmyfT7eZ552sl8c/JNcpK8OTnnTQwzMwEAAAAAqBBrtwEAAAAAOBsECwAAAADQBMECAAAAADRBsAAAAAAATRAsAAAAAEATBAsAAAAA0ATBAgAAAAA0ibfbQDhccMEFVFZWRo0aNbLbCgAAAOAqjhw5QomJifTTTz8F1bo6WCgtLaXKykq7bQAAAACuo6KigvTmZXR1sNC4cWMiIlq5cqXNTgAAAAB3MXjwYN1a9FkAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJrEW1n4oEGD6ODBg9Wm33jjjfTUU09Vm75o0SL65z//6TctMTGRfv75Z8s8AgAAAEAbS4OFBQsWUGVlpe/7nj176NZbb6Vhw4apzlO3bl1avny573tMTIyVFgEAAAAQBEuDhQYNGvh9nzVrFqWlpVHv3r1V54mJiaFGjRpZaQsAAAAABohYn4WysjL66quvaOTIkZqtBUVFRTRw4EAaMGAA3X333bRnz55IWQQAAACAAhELFr799ls6ffo0XXvttaqali1b0pQpU+itt96il156iZiZxowZQwUFBZGyCQAAAIAALH0MIWfhwoXUv39/atKkiaqme/fu1L17d7/vV1xxBc2fP58eeOCBCLgEAAAAQCARaVk4ePAgrV+/nq6//npD8yUkJFD79u3J4/FY5AwAAAAAwYhIsLBo0SJq2LAhZWdnG5qvsrKScnNz0eERAAAAsBHLH0N4vV5atGgRjRgxguLj/Rc3YcIEatKkCT388MNERPTGG29Qt27dKD09nU6dOkXvvvsuHTp0iEaNGmW1TQAAAACoYHmwsH79ejp06BCNHDmy2m+HDx+m2FipcePUqVM0adIkOnLkCNWrV486duxI8+fPp9atW1ttEwAAAAAqxDAz220iVAYPHkxERCtXrrTZCQAAAOAujFxD8W4IAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAJGjuNhavdllFhfr1xvRGvFhlV8neLBy+x4/rl/rFM/Ruj/cts1C0dcAECyAyDB7NlGXLkT5+fr0+fmCfvZsezzMnk3UoQNR+/bBPRjREulfN6v8OsGDldv3hReImjYV/rrFc7TuD7dtMyPl1jTYxQwaNIgHDRpktw0QjKIi5tatmYmYMzOZPR5tvccj6IiE+YqKIuuhqIg5I0PQEjGnp6t7MKJl1r9uVvl1ggcrt29uLnN8vKCNjxe+O91ztO4Pt20zI+VGCUauoWhZANaTlES0ahVRZibRvn1E2dnqUX5+vvD7vn2CftUqYf5Iejh6lMjrlb4zC9PC1RpZN6v8OsGDldt32DCiigqi+Hjh77BhzvccrfvDbdvMqnNPtGBx4GIpaFlwGfKoXSnKD/a7RR5eeIE5MZH5hRcCfk9LEz56/BrRGlk3I9vMbR6s3L45Oe7zHK37w6HbbMUK5vbtmVesCKNcl2PkGopgAUQWtYMykgerbFmVGZlcJ6mSiZjrJFVyZUaAByN+rVq3aPZg5bo5wQf2hyO3mbdlJvfqUsJEzL26lLC3ZRjluhgj19AYZma7WzdCZfDgwUREtHLlSpudAEMENvfNm0c0bpz0ffVqotTUiHh4et9YmkyTfZOfpceo+PzmtCLjdqLEWsLEslKiXbuISkuJatUiatWK6LffaHjpQpqU+ZGf30v7ldDpn36tpvV9b9+eKLEWPfMM0ZAhkp1ffiEaP17Dr8zDiow76bwP3/Zts3eTH6VZLSYH9duhdAvNzXzWz++ddxJt2xZ8c40fT3TXFdJ+K85oT9l1NiquW6CHd5o/Q90+fczn93/NbqSJTd4L6rd26Z/0XeZ4P7/PP3qSvnj9QNDtm51d1b9RVteuqfMNFWRcGNTzxNLJdG3mdl+9PLCvkkbX+kJ1/eQ+vii9nJpl1vF5nj+f6NUXlLVyD+ltatGnn/pv8wcfJFq/PnidGF36Pj2U+aXfcdS31ibytu8Y1O8rpffQxZkFPr8//EB0713B/VJiLVq/nijukLR9pzecTB82eTDo9u1VupbeyHzVz+9Ndb+gva2HBfV7b+k0Gpu5wef3yBGi4cOD118iog8+IGpTW/D7332taRj91/fbchpKQzP3Rubc4yAMXUMtDVssBi0LLkZ+lyB+IhzVl+zxcAxVMJG3yoKX68Sc5VFXnvWzpfa5pe5n1fyef37w+YiYP/7Y38vGjfrmI2I+Tuf7bbN/PXxS13y9am2t5vfii/Ut84knqmao2m9nqI5uv9/TxX5+F75zRNd8dWLOVvN7xx36lnnttbKZqjyn0QFd876T/JjfhF0tLtW9rgdS+/l5fu01ffO1bVu9fl5+ub557ztvbrXjKDbWq2veJU3H+/ldsUJ/Payo8N++/6AXdc03OGltNb+d25XqmvfFBs/7+T14UL/fbduEebx5Hu6VuLXq2GeOoUrulbiVvXk1p0VBBB0cgfNJTRXuLOTMmxfRqH7YHanEFEdEMVVTYqiI69CvnjoUE0MKH6YY8kqfAZdU86uqjWG/spRQXmZ1D37Mm0dUr54+v60yNfwG/xCR335TW7egfpOT9fmtXSvk7eu3jas8xxDr8lytiefFFzXXT+6DXnvNtDqhe98MHeI/47x5FBMTo8/vww+H79fg9o3p2LG638REfX7HjAm5Dot8syuVfizrWnXsEzHF0o9lXembXTWnRSEkLA5cLAUtCy7G5paF779nlloU5B+v0HehMgy/Vq1bNHuwct2c4AP7w1qtTrxeoY9CHJX7FRtLFULfBW/IRbsSdHAEzkQcsxzYUWndOuWOS3rHOBsZC11UxKf+KOL69So0myyffeSEUG5RkT6/RrShrJvZHkStXg9WbAcrt++xY9Zut2iuE27RhrDNls8r1Dzul88rNFauy0GwAJzHrFlCkhO1IW2BJ4ecHEE/a5a+cvXcbcyaxZyRwbfU/ljzhEHk5Tp0hitTWjA3b86cnKztNzmZuYVObSjrpmebGfEgatPT9XnIyDB/O1ilzcxknjiROSFB+GvFdovmOuEWbQjbzLshh3slbuVYUr5RiKUKoe/CBp3lRgEIFoCzkGdRE7PrKTUnyk8Gos7kjHaf0zVBAgXpxPEn1ZUmtGyp7FeedS6YNtR1C7bNjHgI1BrJqmfWdrBKK27fwI8V2y2a64TTtSFus5K4OtyEDmse903pEJfE1QlebpSAYAE4j5wc/zS8OTnh6UQC7zbUAgaPhw836crJJDVD3nF1AX/wAVf/PJ3LP1Av6QwSF6ftNzZWv9bIuhnZZkY8xMVJ2ubNNbcZp6RIWjM9WKUVt5cbPdtZJ9ymDXGbeeIyeNO/t/OmTcybNjG//LK0yGsuOcb5cen6y40CECwAZ2FVy4LSfCrletMz+HJa4jsxjKBFQiIWF9wRmXoXKXqQBwx61s0Nd7JoWXDWurlgmxUVMdepakhoGHuMyynO2LnH5SBYAM7Dqj4LImoBQ9X0U1SXL034HxMxN21YykfSe2p7cNCzVtOfT4se0tP1rVtGhnuekaPPgrPWzQXb7PorpLwqq5rdaPzc42IQLABnYtVoCJEg5VYe8PDr00p52TKdHhzUi9syD3rLtcIDRkM4q064TWvSNvu48X1ci4r56jrfcM6XBcbKdTkIFoDzUWo2lp8EIlGuE7TRvG5O0DrFB+qEtVojBJRbTLX4dEan8Mt1IQgWgDtYt87/RLBunenlllG8drlGPFilNYIT/LpN6xQfqBPWao1gVbkuA8ECcD4RuGtYR305g/YJzyGj8I7I8Xd7TtA6xQfqhLVaI1hVrgtBsACcidUZHGXlnsrozJkpxUwkvChmRdOxznnWGs3ZE52gZXZfnwXUCUdss5I9Hl6xgtl7Fn0WAkGwACKD1RkcA8q9bcxp301D31qbhCFRTujFHc3ZE52gzcx032gI1AlHbLOpNIHPiznFRMxbU4djNEQACBaA9UQig6Os3M9nS69ArluXee+ag84YHx7N2ROdoBXrTuBHra45wTPqhGO22ZsNJ/kkT9LTyLMQAIIFEBmsyuAYoD+8ZJPv5oGIec4cmc4JmeeiOXuiE7RuzOCIOuGIbXZw8SafrBP9jAyOASBYANYToZYFLxFfkbTKd8CPGCG8ktYxd0TiukVj9kQnaN3YsoA64ZxtlpnJfWmdT7477VK0LMhAsAAiQwT6LLyd/LjvQG/SqIILCwPKdcKzVnHdojF7ohO0buyzgDrhmG02rcFzvnPI89du1D73RAEIFoAzsXA0xO7dzHWSKn0H+pKm453bi1tcN73l6t1mTlg3J2iZ3TcaAnXCEdvst+8P+s4hF1wQ2mnOTTgmWJg+fTpnZWX5fYYOHao5z9KlS3no0KHcqVMnHj58OK9evVpVi2DBxSg1G8tPAgb5/HPphTB3nzdPX7lGPEDrTq1TfEDrGm23btLkvLzqs0cTjgoWrrzySi4sLPR9jokRvwKbNm3i9u3b8+zZs3nv3r386quvcseOHXn37t2KegQLLsfkLGq7dzPfeCPz2ZUb9JfrhMxz0FqrdYoPaF2hfeYZafJrr6kXEQ04Kli4+uqrdevvv/9+vvPOO/2mjRo1iidNmqSoR7DgYkxuWQipXGijX+sUH9C6RvvLL9LkSy6pPns04ahgoWvXrnzxxRfzoEGD+KGHHuKDBw+q6gcMGMBz5871m/b666/zVVddpahHsOAyQnnWqkFlpaxcu5+fQuscLbP7+ixA6xit18vcNquSMzKYH3mkakRVlOKYYGH16tW8dOlS3rVrF69Zs4ZHjx7N2dnZfPr0aUV9x44defHixX7TPvjgA+7bt6+iHsGCizAyGkJndrY77mD++8AdfDatnXt6cUNrrTYz032jIaB1jjYzkzknhwsy+rD3HWRwlGNpsBDIn3/+yT169OBPP/1U8XcEC1GKkTwLOrOzffGFJOlFG7mSYoQvTh4fDq21WvGEH/jRW9ecvn7QRqb+6MnxEiU4NlhgZr7uuut42rRpir/hMUQUYySDY5DsbIcPSzcIRMyz6TZ95Toh8xy0yOAIrXO1RrLHRgGODRbOnDnDvXr14v/85z+Kv99///181113+U0bPXo0Oji6nVAyOKpkZ/N6ma+4QvrpavqCvXrKdcKdC7RoWYDWuVqVlgWPh3nnTo5KHBMsTJ06lTdu3Mj5+fm8adMmvuWWW7hPnz6+4ZOPPPKIXyvDpk2buEOHDvzuu+/y3r17efr06Rg6GS2EksFRITvb229Lx3xjKuA/GrRDRjto0WcBWlP7LHDr1nzoxXncu7cwecQI1TObq3FMsPDAAw/wxRdfzB07duRLLrmEH3jgAc6TZbkYO3YsP/roo37zLF26lIcMGcIdO3bkK6+8EkmZoolQMjjKtLtbDOI6tSt8wcLXTccLvyOjHbSBdQejIaANp/4UFXFFBXOTJsJPtWszq/TLdzWOCRasBsGCi5EftOJHfrAGaMtaZnEv2uiT3nXuh6paI+VCG+Vap/iA1pXa//s/SfLZZ9WLcjsIFoA7MJBx7clbPT5ZG9rNZ77dYEq50NYArVN8QOs67TffSJK//EW9OLeCYAE4HwPRfcV+Dw+ovYGJmOOonDdSL0ffjUDrIK1TfEDrSm1ZGXP9+oLk3HOZS0qqF+dmECwAZxJG34IKiuUp9V/kKXcdcOZzTmido2VGnwVow6s/svwKN98sxRRffx2502UkQLAAnEcoGRzd0oMaWudoMzMxGgLa8OpP1WgIMXvsl19KwcKttyqc21wMggXgLELM4OgVj1Anj82G1jla8YQf+NFR11yxftBGpv4E5FkoKmI+5xxhUoMGzOXlHDUgWADOw0gGx9hYLqDGfBGt5fWxFzs/6xu0ztEigyO04WhVzlE33CAVsXKl8uxuBMECcBYGWxa8RDycvmIi5liq4I8b3evcuxFonaNFywK04WhVWhaYmefPFyZfdhnzmjUcNSBYAM7DQJ+Fd+o+6DuuG8f8wX9QI2c/54TWOVr0WYA23PoT0GeBWYgZxD6z0QSCBeBMdIyGyE0dxHWSKn3BwuL3ClW1vgPeCT2ooXWOlhmjIaANr/7IRkNEMwgWgPORH7RVn/KWbbhP9xLfJN87xRS0fgd2kHKhrcFap/iA1p3aKAfBAnAHAVnUnr4tz/e1TRvmM2fUtU7O+gatw7RO8QGtO7UBlJcL/Ra8Xt2zOBYEC8D5BET3OdSb46iciYQOyxs3qmsdfzcCrXO0TvEBrTu1Abz2mtTtYfPmoHLHg2ABOJuA54anV2zg1vH7fcft0w+eVNWqPmOEFlqlE74TfEDrTq0Cb78txRdPPKEpdQUIFoBzCTxYPR5esYI5MdHLRMx9aAOXt2wj6BS0amVAC21Y06GFVkfAcOgQc0yMIO3QQVXmGhAsAGeicVBu387ct2cJ56YOEn5PSxM+eu4YoYU2sG6pDdF1qmdonaMNEjBcfLHUurBrl6rMFSBYAM5Dx8Ho9VbpxINaPMjVDlxooQ3UinVMK/mX0zxD6xytjoDhlVek4p57Trkot2DkGhpLAFhNcTHRoEFE+/YRZWYSrV5NlJpaTRYTQ0TJyUSxsf4Tk5OVy4UWWjmpqUTLlxPFxxNVVAh/ly9XrGuO8Qytc7SpqcK5KTNTOFcNGiScuwK49lrp/0WLlIuKSiwOXCwFLQsuQszgWBWtz57NPGGCyvvhZ80S0rSmp/tlUVMtF1po5UydKmRwnDo1uNYJnqF1jpZZOEcFZHAMpEcPqXVh//7gRToVI9fQGGZmuwOWUBk8eDAREa1cudJmJ0AXxcVESUm0Zw9Rt25ERUXC33XriOrUUdASESUl6SsXWmjlHD9O1KCBPq0TPEPrHK2o19A+9xzRE08I/7/yCtGDD+or1mkYuYbiMQSIHElJVFFBNG6cECgQEfXurRAoVGl1H9jQQhuI3kDBSh/QulMr6jUYOVL6v6Y8ioi32wCoWTz3HNHGjcL/rVsTvfyyvX4AAMAo7doRtW9PdPKk0Drq9fp3jYhGECyAiLFxI9Ezzwj/x8URffABUd269noCAIBQ+OYbopSU6A8SRBAsgIhw9qzw+KGyUvj+xBNEffrY6wkAAEKlRQu7HUSWGhITAbt5+GGiPXuE/3v3Jnr8cXv9AAAA0A+CBWA5X39N9M47wv916giPHxIS7PUEAABmceQIUW6u3S6sBcEC0EYhKYmmVkE/c6b0/yuvELVpo641Uq4qx4/r18IHfNjhw42e4aMahbtPUHY2UdOmRP/4h87y9fowglF9CCBYAOrMnk3UpQtRfr4+bYcOQhfh2bP9fvr8c+GxwzXXEN15p7bWSLmKvPCCcOS+8EJYnuEDPizz4UbP8KHoI7lTU8rdfJq8XqHD4+nTGvr8fOF8qseH3vOukXLDxdL0UBaDDI4WUlQkZDHT8WIVLioSMqSJKc3S04VpAVRW6tcaKddHbq70ToD4eOF7mJ7hAz5M9eFGz/Ch6eNvMW/5ZvvkExW9/L0TrVtr+9B73jVSrgp4NwQIn6QkolWrpDzp2dnqke7Ro8JAYxFmYVoAsbH6tUbKJSLB27Bh0jsBKiqE72F6hg/4MNWHGz3Dh6aP63iB7yfFBE35+cL5U3w3zqpV6kmfjJx3jZRrBobCEIeBloUIEOxNbAqvgS2nOL657gLesuxwUK2RcnVpnfJ6YviAj2jzDB+KPspbtuEGdJSJmOueU8nFxRr6YC0FeucLtdwA8IpqYC5qFVNpusfDk+u/wkTMCVTK7750TFNrpFzd2hA8wwd8RNyHGz3Dh+L0W+t+ykJzBPNX7xZq6/Vi1EcIIFgA5hNYQdetU6ywGzcyx8V5mYg5jso5J+VaVa2Rcg1rrSwbPuCjJnuGj2o+Fr9X6AsWbqn7WXDfejG6PQyCYAFYg7ziih9ZhT1zhrlNG+mnp85/TVVrpNyQtVaWDR/wUZM9w4cfxcXCIwgi5vp0jMso3rQLuuHtYQAEC8A61q3zr7Tr1vl+uusuaXLv3sxl361X1RopNyytlWXDB3zUZM/w4ceYMZJsBQ3W51svRreHThAsAGvQiHAXL5Ym1anDnPvdIdujfUvLhg/4qMme4aMan751hGOokvvRGv6WBgX3rRe0LIQPgoUIII7b1Xh29kd6L26cXOGrxzOnHAv+nK2oSPjoeSZnRMvMfEylUyV8wIdTfLjRM3xo+ijKaM+HqYm23mAeBN0+QgTBAjCHWbOERB9qw5Q8Hva2zOSr6QtfoDC89jfsbZisqPWVkZzM3KIFc/Pmwv9maTMzmSdOZE5IEP6qeIYP+LDVhxs9w0d4PjIzhfNo69bCeVUPgfOrlRtGwIBgAYSPPJOYmDFNoWJuWXaY46mMiZgb0R9cQI0FbcuW1Suxx+OfJU38mKEVD57Aj9LBBB/wYZcPN3qGj/B8iHrxPKon02KwgMCkgAHBAjCHnBz/1Ko5OYqyH977mbPoV/6SrhK0cXGqWs7JYY6NlQ4wM7WiV/Gj4Rk+4MM2H270DB+GfFRWMv/4I7PXG6DX8iuiNxAwIWBAsADCR2fLghiVl1KC+RF8TbjrgI+a5cONnuHDkI+ZM4WnGUTMm5YcNtay4OB3Q5Chkh0GggWL0dFnAc934QM+aoBn+NDtY+aUY7744bHz35T0evssiOddvS0FHo+xvhAyECwA81AYDbExZQRPvi2Py1u28T+g0HMcPuAjej3Dhy4fBWm9OCZGyGLbjnb66/Xe+RsdNWFUXwWCBWANHg+fyejIbWg3EzFfSOu5IK2XcgSs1FwnP2gioYUP+HCDDzd6hg9N/SX0nW+WnSsPKesdAIIFYBn/N+Kw7yDoRRuFLI1qOCCrGnzAhyt8uNEzfKjqX6X7fbM8+6z2LHaCYAFYwtdzpZel1KEzvJva4C4MPuCjJnqGD019HqX6ZuneqVRZ7wAQLADTKdzyOzeJk4KFt/+xV/kZHrP+531WauEDPtzgw42e4UOXjwtqbfOdL/etPchOBMECMI+iIvbmefiaOt/4Kv6Vg4qE8cOBB0hurvBROqgipfV4hE5HStPhAz6c4sONnuHDkI8pE074zpkvN3hWmB5iR0SrQLAAzGHWLOaMDJ7T4BFfpU9uUMGHD8s08gMlLk746Im+rdI2aCCMZ65fHz7gw7k+3OgZPgz5+PVX9p03L6K1zGlpQj6HEIY4WgWCBRA+RUXMGRm8lzL5HDrtq/RffKGgzc2VDjDxgMvNVS7XKq3HI2VCET+pqdUPdviADzt9uNEzfITso2NH5gbnV/D4c+ZzBcUK2owMx7QwGLmGxhMAajDTc/Q4naW6RER0+5gzdM01dW02BQAA7uDrr4la8CGKz55AdNYrTGS211SoWBy4WApaFixm1iwuTsviB8+dzW1oN5/O6OSOJtu4OOEvfMCHU3240TN8hOcjLY05PR2PIewAwUIEqMpmdiajY/UDQn4gGO0cVMM7P8FHDffhRs/wEZ4Pj8cxjx9EHBMszJw5k6+77jru1q0bX3jhhXz33Xfzb7/9pjnPwoULOSsry+/TqVMnRS2ChQgSWPExzAw+4KNmeoaPkHycOMEc5PIXcRwTLIwfP54XLlzIubm5vGvXLr7jjjs4Ozubz549qzrPwoULuUePHlxYWOj7HDlyRFGLYMEaXn6ZeedOhR/kB4T4UTpwnKKFD/hwgw83eoYP3fqj237nYcOEd04NH648q104JlgI5NixY5yVlcU//PCDqmbhwoXcs2dPXeUhWDCfJUuEOl67NvM77ygInJBW1eWpYOEDPqLCM3zo0ldWMqekCF8TE5n//FO7iEji2GDhwIEDnJWVxbt371bVLFy4kNu3b8/Z2dncv39//r//+z/OVRnCgmDBXAoLmZs0ker5W28FCJwQwUfZXQd8wIcrPcOHIf3f/y5N+vhj5dntwJHBQmVlJd955508ZswYTd3mzZv5888/5507d/LGjRv5rrvu4h49evBhv0xAAggWzMPrZR4xQqrQV1whTFN6RbXqMzy8zhc+4CN6PcNHyD5Wzf/Dd24dNSqkU7QlODJYePLJJ3ngwIGKF30tysrK+NJLL+VXX3212m8IFszj3XelQCE5mYUsjbNmMbduzZyTo3ygyA+Q5GQhUUnz5sL/dmgzM5knThQeDk6caJ9n+ICPaPMMH2H5KG/ZhpPPL2Mi5nNqlTlmUITjgoXJkydz//792aPWfBOEe++9lx988MFq0xEsmMPevcx160rBwuefsxBdt24tTIiPr34giHg8QkYycWbx07JlZLXiwRr4ibRn+ICPaPMMH+H5qNLfFvOuT/rFJyXsBBwTLHi9Xp48eTL369eP9+/fH1IZFRUVPHToUJ4yZUq13xAshE95OXPfvlJ9Hz9e9mNOjhQoxMcL35XIyWGOjZUKiYuzRyt6FT92eYYP+Ig2z/ARno/4eF5Cl/vk48YpSyONY4KFp556inv27MkbN270GwpZXFzs0zzyyCM8bdo03/cZM2bw999/zx6Ph3fs2MEPPvggd+7cmffs2VOtfAQL4fPMM/6B8alTVT+gZQE+4AN1Gj7C91GlL4mrw+fRSSZiPv98L5eWsu04JlgITK4kfhYuXOjTjB07lh999FHf9+eee46zs7O5Y8eOfNFFF/Edd9zBv/zyi2L5CBbC49AhYSgPkRBUr18fIECfBfiAD9Rp+AjfR2Ymc04O31j3S19s8d//su04JliwGgQL4bNmjZCu/IknVAQYDQEf8AHP8BG+D2Ze8FEJEzH37Mm8dKnmqTkiIFgAhvjzT+ayMh1CpSY4+YHjNC18wIcbfLjRM3yEpC8qYj5wQHl2O0CwAKzDCZnSoii7G3zAh2s9w0f4eptBsABUOXKEefp05srKEGZ2QgQfxXcd8FFDfbjRM3yEr3cACBaAIl4v87XXCnU4O5v5998NzKz3GZ5TtPABH27w4UbP8BGejyq8Xubt26sy5doEggWgyHvvSQFvw4bCaIigKHX2ESt+4HSr3htfE951Dx81y4cbPcNHeD48Hl+H8ffeY27VSvjpxx/ZNhAsgGr89pt/lsZFi3TMNGuWMP44LU1fRB0XJ3zs1DZoIOSFqF8fPuDDuT7c6Bk+wvORliacT2fN4lmzpHPxxIlsGwgWgB8VFcwXXyxVTr8sjWoUFfknKklLq34giOTmSgeYeMCpvCnUMq3HI4yRFrVEzKmpkfcMH/ARbZ7hIzwfHo90w0XEnJHBfxwo8iWMbNPGvkcRRq6hsQSinhdeIFq3Tvg/M5Potdd0zshslSUAAKiZMFPjxkSXXCJ83bOHaOdOey3pwuLAxVLQshCcn36SsjXHxhocyTNrlpCxyW2PIeLihL/wAR9O9eFGz/ARno+0NOF8OmsWMzO//rrU2DB5MtsCHkMAZmY+e5a5XTupQj7+eAiFoIMjfMAHPMNH+D48HikjbtXP4rm5a1e2BQQLgJmZn3pKqow9e+rM0qhGYMXHMDP4gI+a6Rk+wvMho3dv6Ry9d6+ixFIQLABmZj59mvnOO5mTkph37TKhQPkBIX7UDgQnaOEDPtzgw42e4SN8PTNPnSpJX3xRVWYZCBaAHxp11ThOSKtaE1LBwkfN8uFGz/ARtj43V5JeeKF20VaAYAFYgxMi+Jpw1wEfNcuHGz3DR/j6Kjp3FnLg/OUvzOXlmlLTQbBQg1myREjAZBp4RTV8wAc8w4d5PgLIzWUuLlY47ypoVTGilYFgoYaybx/zuecKUeq775pQ4KxZzK1bM+fkKB8o8gMkOVlIVNK8ufC/HdrMTCEdWkKC8Ncuz/ABH9HmGT7M8SEbOqmKmDlXj1ZcTuvW+rQBIFiogQRmabz11jALLCoSKiCRlKhBfiCIeDz+mR7FT8uWkdWKB2XgJ9Ke4QM+os0zfJjrIz1dvSUgMHOuljbQT+vWhlsYECzUQJ57zr+OnzplQqE5OVKgEB8vfFfTiblLiYREJXZoRa/ixy7P8AEf0eYZPsLzIU8l3bx59aCiijO/5vOBlL6SVivNfmBLh5pOAwQLNYzALI1r15pQKFoW4AM+UKfhI3wfol4eMASUfeIE87XDznJSTBFfQV8LQYLezLkhBgrMCBZqFIFZGh97zMTC0WcBPuADdRo+wveRmSmcR9PTFcuuPODhFnGHmIg5kUr45I589YDApECBGcFCjeLvf5cChR49mEtLTV4ARkPAB3zAM3yE70M8n6po76XXfefyDz9kY+WGCIKFGsLy5VKgULs2886dFi9QXnHFj1qFdYIWPuDDDT7c6Bk+TC97dbMxvq8jR4bhwwAIFmoARUXMKSlS3ZkxI0ILdkKmtJqQ3Q0+apYPN3qGD1PLrlizjhs1Er7WqSM8Yg7Zh04QLNQQli9nbtaMeehQZq83Agt0QgRfQ+864COKfbjRM3xYUvbtY077vi5aFKIPAyBYqEEcPcpcUBCBBRl5duYELXzAhxt8uNEzfFjmY1nTW3zxwNixIfgwCIIFYB5KHXLECho43ar3xteEd93DR83y4UbP8GG5j1JK4PNiTjERc71zK7i0ZVt9PpDuWZuaFixUVDB/9lmEHjkwS2lH9Y73jYuTxhLbpW3QQMgLUb8+fMCHc3240TN8RMTHTTSPiQTZMhoavNy0NOE8jXTP6tS0YOH554W6ceWVEXj0EJh2VCuTWG6uf8KRuDhhWiS1Ho8wRlrUEjGnpkbeM3zAR7R5ho+I+lgUex0TMXegHfxF7AhtH+KNHJFwvrYw3XMsAVeweTPRk08K/y9bRrRnTwQWyhyBhQAAABAZRv+lXdSOfqFOdE3MYv0zWn2+NhSGOIya0rJQVMTcvr0UQP7znxFa8KxZQsYxtz2GiIsT/sIHfDjVhxs9w4czfaSl6X9DZQB4DBFl3HuvFChYkqVRC3RwhA/4gGf4cLYPdHDUpiYEC//9rxQoRCRLoxqBFRTDzOADPmqmZ/iIqI89qQPZq9eHQRAsRAlHjwpJl8RgIWJZGtWQV3Lxo1ZhnaCFD/hwgw83eoYPy8v+6CPmTp2EnzamjNDvwwAIFqIAr1fIDy7WjYhlaQyGE9Kq1sBUsPAR5T7c6Bk+LC17zhzpp0fH5hvzoRMEC1HARx9J9aJBA+aDB+12xM6I4HHXAR/R5sONnuHD8rILC5ljY4WfWsfvFx5F6PFhAAQLUcCJE8w33STUiwULbDSCV1TDB3zAM3zY4iO7b7EvPtjeYpi6D3Rw1CaagwWRDRtsXPisWcytWzPn5ChXUHmlT04WEpU0by78b4c2M5N54kTmhAThr12e4QM+os0zfNjiY8Y5j/qChacfPKnsIydHOE9j6KQ6NSFYsI2iIqECEgkpTQMrt4jH45/pUfy0bBlZrXjwBH4i7Rk+4CPaPMOHbT7yqbnva+fOKj7E83Pr1pZmcCRDJTuMaAsWfvmFOT/fbhcycnKkihgfL3xX04kP14iEZCJ2aEWv4scuz/ABH9HmGT5s89Gno/Taar/Mz3rPzxogWHAhYpbG+vWZP/nEbjeMlgX4gA/UafhwgI8X60/xfW3cmHnFigAfaFkITjQFC/fdJ9WP7t2Zy8vtdsToswAf8IE6DR82+9hDrfzih15dStjbUrY89FkITrQEC998I1WE2rWFxxGOAaMh4AM+4Bk+bPXRkn7zCxiW0xB/HxgNoU00BAvHjjGnpEiV4PXX7XakgbySix95hXWaFj7gww0+3OgZPiJWtjfPw50SfuU4Kmci5jgq516JW9mbp+LDAAgWXILXyzxqlFRXhgxhrqy021UQnJApDdnd4CPafLjRM3xEpOzly/1l4mf5cm0rekCw4BLef1/a8Y7J0qiFEyJ4h0T78AEfNdozfESkbK9X6KMgtiqInzgqF/ouhPkKAAQLLuDAAebzzpN2/mef2e0oCHqeszlJCx/w4QYfbvQMHxHzsXxeoV+QEPhZPq+QwwHBgsPxepn795d2+M032+0oCIGVW6zMStOdoHWKZ/iAj2jzDB8R8+HN83CvxK0cQ5WKgUIsVYTddwHBggv4/HNhhExGBvOff9rtRgO1yq30e1qa8LFTm5mpPtQTPuDDKT7c6Bk+IuqjJLU1N6HDioGC+GlKh7ikZbvqZekEwYJLOHyYecsWu11oECxQkOvEg0A8KOzQil61kkjBB3zY7cONnuHDFh87mg7m+HgvEzE3bMj8ww/MmzZVfZYc5vy0i4IvXwMj19BYArbRtClRt252u1ChuJho0CCiffuIMjOJVq8mSk1V1iYnE8XKqlJMjDAtktrUVKLly4ni44kqKoS/y5dH3jN8wEe0eYYP23z8VN6VKipiiIjoppuIevUi6tGj6nNFU2qxdr5wft63TzhfFxcrl2sGhkMRB+G2loUffuCwe69GFDGDo56IddYs4ZlKenrwTGJWaZmZp04VsqpNnWqfZ/iAj2jzDB+2+Li8U56vQUJ1FKbHE5EMjjHMzNaFItYyePBgIiJauXKlzU6Cs2IF0ZAhRKNGEc2cSdSggd2OdFJcTJSUpF9LpE9vlZaI6Phx/RsYPuDDDh9u9AwfEfVx7BhR05ZJVFFBlJZGtH+/fwNFtbL1epZh5Bpq+WOIDz/8kAYNGkSdO3emUaNG0fbt2zX1y5Yto2HDhlHnzp3pqquuou+++85qi5Zz/DjRLbcI/3/2GdGnn9pqxxhGKmBSkn69VVoiY5EYfMCHHT7c6Bk+Iupj0TIhUCAiuuEGjUBBLNtiLA0Wli5dSs8//zzdc8899Pnnn1O7du3otttuo2PHjinqN2/eTA8//DBdf/319MUXX9DgwYPpnnvuodzcXCttWgoz0d13Ex06JHy/7DKiO++01xMAAABnM3++9P+YMfb5ELE0WJg7dy7dcMMNNHLkSGrdujVNnjyZateuTQsXLlTUv//++3TJJZfQ7bffTq1ataIHHniAOnToQB988IGVNi3lww+lloT69Ynmzg0SIQIAAKjRFBQIfcqJiFq1Ejo02o1ll62ysjL65Zdf6KKLLpIWFhtLF110EW3ZskVxnq1bt1Lfvn39pvXr14+2bt1qlU1Lycsjuuce6fs77xA1b26fHwAAAM7n5EmioUOFgRZjxggDKOwm3qqCT5w4QZWVldSwYUO/6Q0bNqR9+/YpznP06FFKDhhS0rBhQzp69KhVNi2jspLor38lOnVK+D5unNC5EQAAANCiXTuipUuJjh0j8nrtdiNgWbBQ03n5ZaI1a4T/09OJZsyw1w8AAAB3EXCvbSuWPYaoX78+xcXFVevMeOzYsWqtByLJycnVWhG09E5l61aiJ54Q/o+JIXr/faJ69Wy1BAAAAISMZcFCYmIidezYkTZs2OCb5vV6acOGDdS9e3fFebp160Y5OTl+09avX0/dHJvmUJnmzYmuvFL4/5FHiPr3t9cPAAAAd7BuHVF5ud0uqmNpv/xbb72VPv30U/r888/pt99+o6effpqKi4vpuuuuIyKiCRMm0Msvv+zT//Wvf6Xvv/+e3nvvPfrtt99oxowZtGPHDho7dqyVNk2nUSOiRYuEURD/+pfdbgAAALiB334j6tePKCWF6MUX7Xbjj6V9Fq644go6fvw4TZ8+nY4cOULt27enOXPm+B4rHD58mGJl4wh79OhB06ZNo9dee41eeeUVysjIoDfffJOysrKstKkPgxmyYkqKadQokxNluC2bIrTRrxX1RpPCOKEuEzkjEyC0Em7bH6Jeb7lBdJ98Ivw9epSqmhcS9HmIBIaTSTuIiL0bQsc7Eo4fZ/7jj6ovYeTqDseDn9bunPbQRr+WObS67oS6zOyMdwxAK+G2/cGsv/7rrPNdukgvp9yfPsDca4gCeEW1mRQVCTtZ4zWgXi/zDTcwN27MvPi9Qun1pa1bC/NHwIOfNiNDqnHp6eoeoIU2VC2z/6t69dZ1J9RlZubcXOn1wvHxwvdI+4BWwm37g1l//ddZ53fulBZ9Ya1N5l5DVMArqs0kKYlo1SrpNaDZ2UT5+X4SMUtjYSHRX29PoD/3HRX0q1aZk7NbhwcfR4/6D8xlrmrTghZaE7X5+UI9FF9hrreuO6Eu5+cTDRsmvV64okL4jmPKvrrktv1hpP7rrPPiIwgiotGl75t7DTEDy0KWCBDRV1TLo8iq6HDFCiHwq1NHigg/oVHB75hM9KD6e1qa8IEWWiu1odZ1J9TlzEzmnBz7fUDr3v1hpP5rzOf1MrdtVcZEzDFUyb+n9bXmGhIAHkNYhWxne1tm8gVdSnxBAhHzWHrfukBBwYPfspSmQwttJLShUFSkfxm5ucLHbK3Hw3zsmP0+oHXv/vB49D8m0KjzW5cf9l1H+tfOMVZuGCBYsJKqnb2chvgFCo2ogE9mdIlINFitwq1bp34ChxbaSGiNIO/sFWwZU6dKHc6mTjVP6/H4d6izywe07t0fYv012sFRYd0mnv+W71ry1nPH9JcbJggWLMab5+FeiVuZyOvbwSlxh9mbF+LJMxTkFU78qJ3AoYU2Elo9KHX2UltGbm71DmfyO75QtR6Pcoe6SPuA1r37I7DeGu3gKJvXS8Qt6TcmYo6N9fIfm3/XV64JIFiwmOXL/euO+Hn88QgbWbfO38C6ddBCa69WD0pNsUrL8HikZ8dEwv9macXlixcou3xA6979EViHQ6zzZymJ76dXuVnDUr7skmJj5YYJggUL8XqZe3Up4Tgq96s/RMxxcV4+eTJCRpxwJwkttKEiL1feoUx+ohanqf0frjYz079DnV0+oHXv/jByLASp8xWpGfxH8+7mHWM6QLBgBVVNQcvnFfrVm8DP8nmFfnq95erWKnWSUXqWBi20kdAarcNylO7k1q1TvruzSusUH9C6d39YVecjgJFrKF5RrYfZs4lefJF43gc06bbaFEsNyEtx1WSxVEmTbjtEQ1rvo5hxY4kmTCC6446g5dKqVUSpqcE9TJkipAAtLRXG/2ZmEq1eLcy7erU07rdHD6LatYWqBy20Vmmzs4k++ohorI66DgBwNxYHLpYS6QyOJXF1uAlJQ1yUPk3pEJfE1WEvEW9PvSLsrF4+rbzjDRFzy5bV5/F4quughdYKbeAzZqMdsfAYAtpo2B9hPob4Hw3gQ9Q0vHLDAI8hzCYnx3dS9MRl8KZ/b+dNm7j659/bOT8unb1E/EDMa5wQX8lffqlRbmAzr1rl8HiYU1KkyhofL3hS8xobK2nj4qCF1hqtvPe6mk4JpXqPDnU1W+vG/RFmB8fy79ZxYyrgGKrky2kJe9eGUG6YIFgwE3kLQGBFliPbyR/GjvXVq4QEL3/1lUb5wSpH4B1fMA9OuOuENvq1obYsYOgktNGwPwLrbQhDJ79tJl0nRtJnxss1AQQLZiMm1FBLRRpwwa9Yl8M31f3CVxESEpgXL9YoXy1gkE9PThYqsh4PycnMLVowN28u/A8ttFZoxWbjcN46GVheYCdKJAGqGVo37g+x/oaYlOl2muW7Rnza+J7Qyg0TBAtWIEZ3wSpx1Um14nQR33ijFIwmJjIvWaJRvp5yjXhwQk95aKNfKz82jB5PeoLkzEykF64pWjfuD/l52UCdL6UEbhB7nImYzzmH+ezu/NDLDQMEC1aj1jwmVqgqysuZx4zxDxiWLg2/XGihdaTWCGoXBKXfg3X8ClWbmem+FxdFs9aN+8PIsSCbb0nT8b7D6S9/CbPcMECwEAmUOr0oUF7OPHq0JKtVi3nZsvDLhRZaR2r1oPek6PEodzgzQysuX6sPUCR8QOve/WHkwh6gH3fdGd9i/DrBRzhgQLBgNQbvtsrLmUeN8g8Yvv02zHKhhdZpWj2EM2Q4PV27E5leLbNyh7pI+4BWwm37g9n/2DDQwbE418Pnnit8rVePuaQkxHJNAMGCFRjssxC4g8vKmK+/XpLl5YVQrhOeUUMLbZC6HhR5Zy89WrHDmZ5OZHq1zP4d6uzyAa2E2/YHs1CHDXZwXLRIikduuSXMcsMEwYLZGBwNodZDvKyM+f77ZedII+U6ofc7tNDqrOtBMRJgiMGL2VpmoUOd3T6glXDb/hD1BnTyx9Kaj6TRwdE8nJpnQdfYcyPlOmFcPbTQyrVG6joAwMd77zH378/cqJFwA2knCBbMRm+2OqNZ7XJyuDjuHB5Jn/G3sZc5P2MftNDKtaFmcAQA8OnTdjtAsGAuFrYsFGd24GG0lImYk+gsr2p2o3PvJKGFVq5FywIArgfBgtmY1GchkNK35vBVdb71nZeT6Cz/r9lfnPuMGlpozeqzAACwHQQLVhDmaAg1Sk4U8fDh0o1cHTrD3zUb7dze79BCG2JdB6Amk5PDvGYNc2Wl3U4kECxYjfwkKn6UHk3opKSE+corpaLOodO8hvppl2vEA7TQRkILAFDl8suFw6d5c9nQeZtBsBAJTM5qV1LCfMUV/gHD93Sx8zP2QQstAECTo0el7j1pacxer92OBBAsWI1Fd1vFuR4elrTaV2RdOsVrU0Y5904SWmgBAEGZNUs6fB55xG43EggWrMCiPguB5RZTLR6a9J2vYl1Ji6uXq/e5s16/Tnj2Da07tUbqulK914sbk/VEu2exfKMYPTcG09ulC9RrMHiwFCz89JP++YwuxygIFszGotEQauUW53p4yBDm/n2K+XRGJ/9y09P19WjPyEBmSGit1Rqp60r1Xm/rhBvTAEe7Z2Zh/1m17+XnRq1l2KUT0bENDh+W0pe0alX1CMLoMRDKttYBggUzsSmDY1ER85kzAeXGxUnhqZ6x8sgMCa3T8izI672exxlufMFQtHtm9q8DZu97pXOj0jLs0hncBjNmSLvh8ccNbAeDywkFBAtmY2EGR93lxsXxH9SIf6Iewl2fWgXzeJhTUvwDBqdnAoTWndpQMzgGtk5o1eXAoCSY3kmvLo5mz3r2Xzjzmn3ONVtnYBv06yfthu3bDc4fzrbWAYIFM7GpZcGv3IwMLqDG3IF28Hl0kjc2vVr/HYXT71Chdac21JYFpTKCHU/i4w69+rQ06UJpl7YmeA714mVk3xs950ZaF2QbeDzSodShQ8AoCKPHgMmBAjOCBfOJcJ8FpXLvqj3XV+nqJZXwDz8EKTcjwz3PvqF1pzbUPgsiaidDM6Y7QVsTPIeK3mUYOefaoQuyDV5+WQoWJk8OYztYECgwI1iwhgiNhlAr9+zufB44oEIKGOox//hjkHLt7ikPbfRrjdR1JfQuIxS9E7RO8WGl51DRuwyn6zTYuJH5jjuYGzRg/vXXMLeDBSBYsBr5zhU/ZuzUIOWeOcOcnS39dP75smE4ZvmFFtpQtaFidBlOWFd4Nge9y3C6LghBX0UdiW2tAIKFSGBVVrsg5Z45I7wLXfy5fn3mTZtM9gsttKFqQ8XoMpywrvBsDnqX4XRduNiQKRXBgtXY1LIgcvo08yWX+AcMmzeb5BdaaEPVhgru0p2nDUUfCk5vMYjUHT9aFqzFlmBBvlMzM8N/vhRiX4jTp/2H5DRowLxvX4jlOuHZN7Tu1MrrmpX1Xsws6KbtEu2erd73Yq6H3Fxn6jS2wW+/Mc+ZEyQhplX94AyAYMEqAnequBPVpgcjzFEWp04xX3yx8PP48bJXnxop1wm96qF1p1ahTlpS7ydOFDILTpzojvpcEzxHYt/HxwtXXKVhjE7QaWyDp54SJAkJzJ9/HuZ2CHVb6wDBghUECwiMBgwm5W84dYr5xRdlgYLB/A1M5P9x8th+aJ2j1aiTmoRS7wM/Tq7PNcGz1ftenk1S/ChllbRLF2QbeL3M7doJk2NimA8eDHE7hLOtdYJgwWz0BgJGAwYLM0NWxCXqK9cJmQChdac21AyOodR7+cnbLdslmj1bve/nzPH3M2eOs3Qa67F1qzT7gAHKs1t27jcIggUzsSqPt4WZIbemDuf29Av/HNfV+Xc10LpTi5aFmusZLQua22DiRGn2t94KYzugZcE8Ip7BMVigIOLxmJ7BUe9zqx07mBueU8xEzI3ql/GOFkO1y3XC81Jo3alFn4Wa6xl9FhS3gdcrxGREQmPOH3+YdAygz0J42JLB0Wy9yT1iT5xg7tVLimwbJ1fwLy2GaJfrhJ7Y0LpTG8qxEUq9j/aRBW70bPW+d+FoiI0bpXPvZZeZtB3C2dZBQLDgVpSaBOWVRScnTjBfcIFURJNGFbyzxWX6yjXiAVpozcDoMpywrvBsDnqX4XRdFQ89JMnefdeC7WAyCBbcjElZvI4fZ+7ZUxYwNCjlXdRWX7lOyB4HrTu1oYJsiM7ThqIPBadnZtSpq6wUnvgQCU9zjh8Pst6h+jERBAtuxeTo8vhx5h49pKKa0iH+lbK0y3XCXQ207tSGCu7SnacNRR8KTm8xMLAN1qyRJMOHW7QdTAbBgpuw6rlVle7Y9t+5e+LPvvrXrGEp/9pisH+5TnheCq07tUbqZDj1Ptqf/7vRc6j7XkTvMpyuq6KggPmll4RHwPPmWbAdLADBgluwqkdsQLnHqD53S9zhCxjeeu6YVK4TemJD606tkToZTr2P9pEFbvQc6r5nVg5elJaRkSEsQxxe4DSdPCALwOs1cTtoLCdcECy4AQvzLCiVe3Tb79ytG/Mrr8jKdcIYb2jdqTVSJ+Ugz4L7PYe675mFwCIjgzktTXsZYnATFycFPk7SidsgLU1YH6MBk97tEO5ygoBgwS1YmMFRSV9SoqBzQvY4aN2pRQbHmus5lH1fVOQfxKSlKTexFxUJLSJyPy1aVA9I7NIxC77FCz2RsF5GHhHr2Q7hLkcHCBbcQIRbFtTuPpbSMN5LsjsLJ9/NQuscLVoWaq7ncPZ9enpowULz5vou7pHQiduh6iJeSgn82HkzeNO6Yv2PH/Rsh4DlMJEwXzQFC/n5+fzPf/6TBw4cyJ07d+bBgwfz66+/zqWlpZrzjR07lrOysvw+kyZNUtW7OlhgjlifBbVyv6SrOIFKOTXud/6tSV93PCeH1jla9FmouZ7D2ffp6VH1GGJJ45t91/L77zd5O8gfQ6SnR99jiO+++44nTpzI33//PXs8Hv7222+5b9++PHXqVM35xo4dy0888QQXFhb6PqdPn1bVuz5YYLZ8NIRauRUU6zdKIi21kvetPRjcgxN6bUPrHK2ROhlOvY/2kQVu9BzOvg8sXykgcUkHx3E3lvvOo19+acF2kO9Lk7E9WFBi9uzZQU2NHTuWn332Wd1lRkWwIEdeQcSPUsRpUrl/bP6dO3SQJqWlMe9fd1C/ByN+oY1+bagYXYYT1hWew0dvQOJgXXEx87nnCpPq1VPoF2bmdrAARwYLr7zyCl977bWamrFjx3KfPn24d+/efOWVV/K0adO4SCOairpggdm6LF4q5RYUMLdvL01OT2c+sPAn/R6ckGkOWudoQwXZEJ2nDUVvFL0BiUN1ixZJk2+5JQLbwWQcFywcOHCAe/TowZ988ommbv78+bxmzRr+9ddf+csvv+RLLrmE77nnHlV91AULEW5ZEMs9fJi5XTvpp4x4Dx+gtOAenHAHBK1ztKGCu3TnaUPRh4rD0zhr6UaPliYvW2ZgncPxYyKWBQsvvfRStQ6IgZ+9e/f6zVNQUMCXXnopP/bYY0YWxczM69ev56ysLM7Ly1P8PexgwegzIAueGdnVZ8Gv3KIiPryviNu1LvPV05bxeZy38EdFre3PVqF1jtZInVTC6DLctF3QZ8HY/tcKSByoO3OGuU4dYVLDhsxlZRHYDiZjWbBw7Ngx3rt3r+ZHPuKhoKCAhwwZwo888ghXVlYaWwtmPnv2LGdlZfGaNWsUfw8rWBB7ZOvdGR6P+e8Tt3k0BBP59Zg+1KAjt6Vdvrravj1z+T5lrSt660NrrdZInVQisCy1ZUydGlJ9ds3IAjd6Dmff6w38HP6K6o/fOOI7V955p7FNYGg7mBWYKeCIxxBioPDggw9yRUVFSGX89NNPnJWVxbt27VL8PeRgQT7WW0/0Jt+ZRrKV6fVgY54Fv0iWiA+l9uaszDJOTGReskRb6+g8ANBaqzVSJ5VQCxS0lhFCfXZ8zgI3eg5n3xu5kYmPZ54zR3m/O0A3os5/fZtq1UOL9a1/KNsh3KBcA9uDhYKCAr7sssv45ptv5oKCAr+hkHLN0KFDedu2bczMnJeXx2+88Qb//PPPnJ+fz99++y0PHjyYb7rpJtXlhNWyEOxkZVQXChHO4KioU8jydvAg87ff6tMaKRfaKNOGmsExlGNP/ERjNkQ3eg41g6OeG5ncXGU/4h2/3TqPh09mdOFEKmEi4W2+Fa2yjD0ituJGMQRsDxYWLlyo2qdBJD8/n7OysjinqqIdOnSIb7rpJu7duzd36tSJL7vsMn7hhReszbNg5O7G7EDBoS0LWncqZ6iOvXdA0DpHG+pJzGirntIJPJru0t3oOZwLmN5AY84cfy9z5jhKV/r9Rv4q9hq+kT7gf8Y8byxYZrbuRtEgtgcLkcKU0RBqAYGVgYKIw/osaD3XfIEe4Taxe/n3Zhe455k6tNZqQ20eNdpfaOrU6H7+70bPoez7KGpZCCtgQstC5DFt6GTgQRDBpBhOGQ2h1WN6Jt3pO3ay2lTywb0u6q0PrbVaI3VSqX7qJdpHFrjRcyj7MYr6LIQUMIWyHcJZThAQLISCUnOc1YFCpDwYKVdB60m7mDPTpJSmWVnMhw6FXy60UaaNFE5YV6PbxQk+rPRshCgZDRFWsGxkO4S7HA0QLISKDUkxIuYhzCxveXlSqnQi5rZthWROjsg0B61ztJHCCeuKDI7hozcgcZDuyfNf549mHGWN7nTWbQeTQbAQCk64g3Joy4KoPXDAv99Uu9ZlfDitd9jlQhsl2kjhhHVFy4J5ODAzo5ruMDXh2FgvEzF37mx8VTWJpgyOTsO0DI42NgNZ7kFPuQY87N9VzOlplb763J5+4YK0Xs59DguttdpQ6mQoIIOj8zybse8d1GKgRzeD7vHJH388vFUPyY/JIFjQgxM6mFjtIXB+pXKTk6VXsurRtmjB+5pcyGmxHl+d7tCmjAsKlLWu6dkPrXGtVcdFIBaM7rF8u2A0hDou7rPQr/YPvvPe9m8O+6+PURxws4pgIRhOGLpitQe1QEH+u/yZgsGcDL9RS06NO8hEzI0aMf/yi7qWiZydMwBa41qrjotALMwbYul2CfxEk+dw9r2LR0N4PNKm6UA72NsyjIDJCTerjGBBH05IimGVh2CBglyXkuIfMBjICLd3wRbu2ZN5x47gWsdnsIM2chkcjRJmRtKIa5HBURmX51mYNk2aPLn+K6EHTE64Wa0CwUIwnLCzIpHBMdgzr6Ki6ncgBu9UvHkOvVuCFi0Ldm6XwE80eQ5n37s4g2OvXtLkXz/ZGl6w7ISbVUawoA8nNANZncFRK1CQazMyTHsGWtGyNT9FT/GRBln2P4eF1lot+iyobxf0WaiOi1sW9u6VJnXvVBpewOSEm9UqECzoxQEdTCzP4KhXa0Lv6oqWrXkc/YeJmLt2KOWj+Q7o4Q2ttVqjdS1UMBrCeZ5D2fcu7bMwZcIJX7AwtcHU8INlJ9ysMoKF0FBqjlOK9qzEbR4CtPnUnJvFFfhm69aN+ejR8MuF1iXaSOGEdTW6XZzgw0rPRnDhaIguiTt9m2A/pZsTLDvgZhXBQqg4IVOd2zwEaH/9aBM3ayZN6t5dullxRFY6aK3VRgonrCsyOIaP3oDEZt2vlMWTaRKPpfetCZZtCsoRLISCE+6g3OZBRbtr1SFu2lSa1KMH8/Ht+fbfLUFrrTZSOGFd0bJgHg7IzGiKLlyQwdE6ouKtk271EES7a9UhbtJEqvc9E7fzCaoXdrnQOlQbKZywrka3ixN8WOk5HFzSsmB5wISWBWsxJVgIPCjEnaM23Qrc5kGndufKQ9w4ucJX93vV2sonfs4Pu1xoHaaNBEqd9NQ85eZWfy5thtbjMaY14tlI2bm51mwLj0d4ZmjEczjP0fUGJE7XhYuNQTmCBb0EO/FF4sToNg9pacJHp/aXZoO5MUmdHm++2ZxyoXWINhIBgzi8V6+nuDjhY7Y2JUX/UDdx2KJez0bKjo+XkqmZuX4NGghl16+vv55kZISW7llPQCKmoRdfd2uDbgbdw0/Xf5V3rTpUXWdFwBThoBzBgh707gwrd5rbPIgnPvFEoVO7o+lgbtSwgjt3Zi4sNK9caG3WRiJgCEwcpuUpN1e6MIoXysCx9KFoPZ7qmU61tIEJkYJtR71lB+YBSEkxZ1t4PEJgI/ecmqq/nmRkGEv3rCfwE/NCiOuQnBxxnZeIW8buZyIhGabiyK5QA6bAcmwKyo1cQ2OpJlJcTDRoENG+fUSZmUSrVxOlpiprU1OF3zMzBf2gQcL8Nc1DcjJRrKy6xMQI03RoO9baS6u/KadVq4gaNTKvXGht1FpVJ5VgtqZcvZSUEBUWSt+bNCFq0UJZm5pK9M03RHFx0jQt/0bKbtGCqGlT6XthoTC/3ejdP8XFRFOmEB04QOTxEKWlKZ9zkpOJatcW/q+sFP7Wrl29Llqs+5F60X5vBhERDR5M1LBhlU6s+2lpwnocOCCsl5H6n59PlJ2tfe4NPMays4X57ML0UCWCmJLBUW+05vFYl8HRLR7Eu4L0dH1ZJINoT51i/vNP88uFNkJaZmvqpJKn9HT7H0PExwvZEPWs69Spgr5BA3PLnjVLaAXQ+8jC6GOIuDj9ntPS9NcTZqH1IT09eItLUVH1lo7mzau3Xlise5Be9smrZYcObF1JTw8tg6OeFgP5NkcGx9AwLYOjVfpo9CBmhAtT++efzBddJHxOnTKvXGgjrBX1VuOEDo65ucbW1UhnQSNlFxVZ14HT6g6OegM/mx9DVMbGc3PKZyLmhASvlCsmcDsYDZjk28HuG0VGsABcwFVXSUH5xRdXBQwABCPwguWGYaJWlW3l+lm5PVzQwXHNzF9856fhdb41N2CSbwcr9TpAsAAcz5YtUmsnEXO/fsynT9vtCrgC+Qlb/KhdwKzSWunZqnKNerByewSWrxWQ2KC75x5plefRTeYHTA4BwQJwBZs3S6O0iJgvuQQBA9CJU9IhW+XZqnKR7jmorrycuXFjYVLtWpX8Z0YX6wImm0GwAFzDpk3+AUP//sxnztjtCjgatCyEVq7TWhZEHJbG+dtvpUnXXWegPBeCYAG4ip9+Yj7/fOlYHDAAAQMIwAmvqJb7MNuzkbKt2hbM1r+iOhAHtizcdZc06dO3jkQmYLIJBAvAdfz4o3/AkJ2NgAFUIfYcz8kJ/gxazJzYvLnUA94MbWamsHy9PdKNeDZStlXbIjOTeeJE5oQE4a+ZngNx+CuqT2/Zwx9/zDzq8lN8NqOD9QGTjSBYAK7khx+Y69UTjsWrrmIuKbHbEbAd+Zj0YLkFAjMnEgk938PVihcIcfnBxrob8WykbKu2ReCds947cr3bQ46RYCc+XkhwoLSuTtCFEzA5BAQLwLVs3Mg8bhxzaandToBjyMmRTtzx8cJ3NV1srHSxi4szT6tn+aF6NlK2VdtCnkbabM8ieoOdwLTWoj4wXbVdunADJgeBYAEAEB2gZcH6bRHJlgW9gcacOf5eqqVQtFkXasDkMBAsAHcTcOLZv595/HiV85HdmQmNdkrTW6bZ62TV8vWWGQ7os2D9tohEnwWHtywsmfU7/yXmI/6SruISSkTLQgAIFoCzCEiDum+flE5+yBDm4uIArZ3vPNCbslV+gtdTptnrZNXymfVtJzPAaAjrtwWz9aMhHNxnYfQ5i32xxLIHlqHPQgAIFoBzUHjByoYNzHXrSjcEw4ZVBQyBry7WepGLES2z/wlB7Y5B78tglO6mtMo0e52sWj6zvu1kBXqHx1mptdKzVeUa9WCVZweOhjiT0ZHr0BkmYm54fgWXlekoD6Mh3AOChSgk8ATi8fDatcznnCOdry6/nLk41+P/1je1N9iJZRrRKp0QwtHqfb5p1KderVXL17udrAAZHEMrFxkcFXUf02if/K67QijPpSBYAO5G4UL0/ff+AcMVSauE54ppafpfXWxEq/eEEGwepeebZvo0orVi+XacONGyEFq5TmlZCMQBGRxH0CKffNWqEMtzIQgWgPtRuCCtWcN8Tp1K3zF7ZZ2VXLLHo37xUppuRBuGV8Xpeju7GfVp9/IjAfosWL8tmGtkBseTdB4nUgkTMTdtVMEVFSGU51IQLIDoQOGEtbrZGN+zRSIheVNpqbJW9aJmRBuGV10X20i/Ytnu7RQKGA1h/bbIzIxcBkcRs+tsiLr/PJHrO5/ce96/QztWXAqCBRA9KET2/2v2F66TJLUw/Pvf6lrVg9qKOwaz736s0tq9nYyAPAvWb4vA/au37oQ6bFCppUPpAp2RIQQkLVtaqrv8cmmV11FfKRAKVl4UBAwIFkB0ofDMcNUq5qQk5kceYfZ6tbVGyrXCa9jLtkJr93YyAjI4hlau0zI4MktDc4P1iRFbQuLihL/JyZbojh6VViWteTlXtmzlHwiplRclAQOCBeBu5HcpGne2e/Ywe88WKT/DVbo7KjKgDfQRple/pEhm+tSrtWr5erdTqKBlwfptEamWhcChuWqjbYqKhMcnci8tWlRfhgm6Dz6QJv3jH1w9cFIqL3A7uDhgQLAA3Is80VGwZ4ZTp0oJhKZOZc7M5B3UgctaZunSqpYrLttIUiYtr8nJUhNoerr/8+NwfOrVWrV8vdvJrDqBPgvu7rNQVCRlWDMaLDRvri8IMKjzepk3b2aeMIF528bi6kGWVq4R+fZABkdng2AhypDfOakN31O6s6lqXlxLF3HdmNM88vKzQlIVDa1mufJl60nKpOVVfvIRlxt4lxeqTz1aK5evZzuZBUZDWL8tmCOTwTE93TGPIVR1KSn6sphGIli2EAQLwL14PMETAwVqiPgM1eFGsUd8k0aNYilgCNDqLlcrKZERrykp/stOSdG1TkbW39blB9tOZqPn0YjVWis9W1WuUQ9WeXZYB0dFnRho6V0fl2LkGhpLAEQB51ARzWv4INWqxURE9NlnRGPHElVU2GxMiZiYmr38cElNJZo3z3/avHnC9EhpjWJV2Vaun1Wek5KEMlavJsrMJNq3jyg7m2j9euHvvn3C9DVriPr0Ifruu8jrUlMFn3rXpyZgceBiKWhZiDLCeQxR9f+yprdwrVpe3+TR5yzmcopT1Dr+MYSeF+ME01q5/Eg+hpCDloXQynVKy0IoyzBZ58k5yBfU2sYv0j84j1KtWTcXgMcQwL2E2cGRiXhp01s5MUHKwzDmnK+4/LkX0MHRrR0c5eh9lm6l1krPTtgWVm+PQGxI9zxtmvTzv+gJ7fKiGAQLwN3oHY4oahU6fC2hy30pXImYb7yRufwUhk66buikHD3PueVBjBVaKz07YVtY6TnYekSwZaFXV+ncsJvaoGVBBwgWgPMJMYHQ13QFJ8RLLQzjxoWRwMlsr0jKZIxgFyq1xyNmao1eTIx4NlK2VdsiM1N9SKYZ20PE5ldU7/12v68qd0/8OSpfPa0XBAsgegjzmeziJrdxQoKXY2JkaaGNlmu2V7ufuTvlOb5e9F6gPB77X1seimcjZVu1LcQytZI9hepZjpE8EfHxzHPm6Ou7Y0A3JfYx32aZOvGE8jqF+94Ll4BgAbgbk8eRf9Xkdv7PK0elss18Hmv2XZLOdTKstWr58m1gBfKOpMH2TWCGQK2EOka0zP7bwkgGR6MBgN4MjmZuC2ahXshTOIv1JVzPav61ghK5F/Gj5ClEXRfa6pPv36+yTqG+98JlIFgA7sXKDHVKHfzCeR5rxV2SVRkIrVh+pO7A5B1J9WjFDpp6siHq1TILyzeawdHIowUjGRzN3hbMQofWhAThrx6MbA8Rve+VmDPHPwiYM8cU3U5q55NeeGEY/qIEBAvAnVid+142dPA/yQ/xbWNOc2WlilbPnZuZd0lG10mv1qrlR/oOzGjmRCMJdULJymi21ojeqm3BLGVwNFK+Ea3NLQtP0VM+6auPFfrr0LKgiWXBwsCBAzkrK8vv884772jOU1JSwk8//TT37t2bu3Xrxn//+9/5yJEjqnoEC1GIVW/VkwUK/z7/fo6JEXIxjB/PoQcMZt8lWfXWRKuWX4PuwIBJ2NhnwRsXz22bnmAi5hiq5N/T+trbYuYAHBMsvPHGG1xYWOj7nD17VnOeJ598kgcMGMDr16/nn3/+mW+44QYePXq0qh7BQpQRoTcMLki+i+PipMRNt92mETDoScqEloWovwMDJmLTaIh9K/dxQoLwtX/tHPXyItEXxyE4JliYO3eubv2pU6e4Y8eOvGzZMt+0vXv3clZWFm/ZskVxHgQLUYjVbxisen772Wf+CQ1vv10hYNCblAl9FqL+DgxYhLw+iR+1GwQTdMePM7/7LvPi9wr1lRflOCZYuOiii7h37958zTXX8OzZs7m8vFxVv379es7KyuI///zTb3p2drZq0IFgIUqx+g2DVXz6qX/AcOedAQGDkaRMGA0BQGjYkMHRkC6KcUSw8N5773FOTg7v2rWLP/roI77gggt4ypQpqvqvvvqKO3bsWG36yJEj+cUXX1ScB8FCDUDvHYVRbRXz5/s/pr/rroCAwQqvVq2T3csHwCgRblkwrItyLAsWXnrppWqdFgM/e/fuVZz3s88+4w4dOnBpaani7wgWgCoWZxv8+GP/gOHuuwMyPVrhFRkcQU3Fhj4L3pYGW9fkPqMYy4KFY8eO8d69ezU/asFAbm4uZ2Vl8W+//ab4Ox5DAEUidBf84YdSwFC3bvCcND7c8m4IuV89ZRrZpjXgpApMwobREJUtW3EX2sq3x8zhVf9YEry8GtQXxxGPIQL58ssvuV27dnzy5EnF38UOjsuXL/dN++2339DBsSaj95m5Ua0KH3zAXK8e89q1Ov255a2T4vp7PPqWr/etn6KmBpxUgQnYlGdhzYce38/D6Sv18mrgKB/bg4XNmzfz3LlzedeuXezxePjLL7/kCy+8kCdMmODTFBQU8NChQ3nbtm2+aU8++SRnZ2fzhg0b+Oeff+bRo0dj6GRNRS3fgdJ0I9og6M5JIz/xqb2sJ3A4orw3pXxIotJJSq7XuqPSq01LY05J0b98rXWSlylqjJ5Ujb7RUysFcqQ0anq9GE26pNeHEb92lSnqguXnEMvSmxvk7beD6u65R/p5Ht2kXV4Nyx9ie7CwY8cOHjVqFPfs2ZM7d+7Ml19+Oc+cOdPvEUV+fj5nZWVxjmxniEmZevXqxV27duV77rmHCwsLlRbBzAgWopZgF3m1C1YwbQgdmLxe5q++UunD4PEEf1mPx+N/kSYSvivp5GUFKzMcrZpOKQjRU6bWS4qU0JOyOLDVRqn1IpKawPU32pqiN02zvJlez9BdPemc7S6TWapfjRurtyyIZa5era9lYfXq6sdAgK68nLlxw3ImYq5NRfwnnYuWBRm2BwuRAsFCFKL34m7kghViwOD1Mk+YIMx2//0KAUOowULz5pENFuTL17pbMhJYhBos6HkZklqrjfzEHUlN4LqLdSmUFyhp1T+lZnqtpGB6XhRld5nM/tusUSPmVq2q91nIza1eplZfBHkLAJHQwqCgW/HhHz7JdTEL9fWBQJ8FRRAsAOdg5Vv1QjjJ//ADc0yMtIgHH5QFDJF4DKGnY1cwrXz5Ri8WwZYf6mMIPcGb3kAsUhq9vsNZZ2b9zeBGAja7ywxcb7XREIsW+Ze5erWgCxy9oFe3bh3ffu583+p8+vohVR1GQwQHwQJwFla+VS+E5uN33/W/0X74YVnA4JYOjvLlG2mG1ttpMpQ7MD2PmpwULIT5OEv3OgcGgmY+hrOzTK2gUNToDVLlusAgRaYrpQSuT8eYiPmcOpXs97YBpfJC3a8uBsECcDdWvlUvhLuFwL5WjzwS0MIgonUCsnvoZKgd3PScVEO9A1O7mOi5cNmlCfeComedxWZwPTqlQNWpZaoRmMdj0SLleRctqq7TKG8JXe6T/uUvOpZbA/OHIFgAwGRmzfI/rzz6qEqnR7uTIllxArTypKrVIqN24bJbE4l1NqJzU5la20IejAYGIfJHD1otAbLyxtF/fNIv52i8jhotC7q0CBYA0Mk77/ifWyZODAgY9J6AjJyorNLqJRInVT3LcJomEuts1ItbylSaRynIkAcM4ic+Xr3lQVaet2Umd8g4y0TM9egEl7RsZ05wE2UgWADAImbO9D93vTSlTPgh2AnI7hdJyR9bBMPIi7z0lhkMPa0XTtOEi52tUHaWqVRflS7kGRnVm/TERw+BupwcodOwrLzKSubvFxTwvEYPauqqlSc/tqIcBAsAWMhbbwnnlPSGp3hferY7XlEt3qmZ+dptvWUGw2mtBmhZsK5MsQNtsA6UYn0O/ATWRVEnjjRKTtYuL5hOXIe0NClYiWIQLABgMf+eVcr70wdIgYDaiU9v6lqPx3/YoviRD7EMVSsPWIi0h07qScdrpMxgaLVeZGQE70cg18jvGK3QhNuaYufrzO0uk1lIjyqvt0qjTsRlN2vmX7cbN5bmlfts3txf16JF9X1TVCRMD6YT9418ZExGRlS3MCBYACASBIw1r1inMtZcb+ranBz/11/GxamPXzeq1ZvCVq/WjLS4Wk3R4p1gSor6CAW5RrxoyO8YzdSE25pixQuUjLQs2VmmuM0yM5kbNlQPFuQZHOU5SYgEP4FBiYqu/JSOYKF5c33BQrD8LS4HwQIAVhNwF15KCXx1nW948kMnecUK5vbtmVesYLQsqKEWKIg+Atcv8MKipAk8uZulCXedzX6BkpH9b3eZgdusUSNh2wYGgPJtJA8AAh8byMtU0V13+Vm+5BLmN99kLinh0B9D6M3f4mIQLAAQCaruhCrX5/DVdb7xnbdSmgi56Ht1KWFvSwvv7NzaZ0ErUBAJTOer9D6NnBz/C0ZsbPUWDrM0gZ6MtqbondeKVii7ywxcb7UOjoH7ISVFuUOihu4knceJVMJEzE2bMlfsly0HHRyrgWABgEhRdUJ59cnj1W7KiJiX0xDrnhkb1cr86lknU8uUl23k3RDygMGKVgO0LFhbptY2k//esmX1/gdia1LgHX/gu1bEfeXx8L8b/cM3+d7rD+kLeDF0UpcWwQIAJvHypOoBwzkxZ/i1p47zoUMyod6e41Zq9WJFmUbeOpmTIzVZW9UfAX0WItNnQWmbKdWv5s2rP6YI7Esg6gIeFVw+sMj38zrqq15frajXLgTBAgA2sHx59Rsz8RMTwzxoEHNZVVqGGp/B0WgLh9UjHTAawtoytbaZUv1SWnZgi8K6dX5lHj0qxS1pdIArKUa7viLdM4IFACKN1yv0UYijctWAITu7Siw7EZaRRrN0gDZqWhas8hFJTSTWx6gXt5Spdz6l3zTKl+dv+ge9aI3fKAPBAgARZvm8QtVzGpHQ6XHmTPY7SVW0bM0ZzYr5qjor+EP6C5/K6Kx+0Qr2bNWK57BOe7aLDI7BvbilTGZ99SuwzJQU1Xo46KJin+yn5ldb8y6LKAPBAgARxJvn4V6JWzmWKhQDhViq4AsSt3LZ9/7Pq1d/+oefrjYV8fXnLOHP3j7CRbn5+k50gXdgwbR6saLMcEDLQnS1LOipX0qdGZU6PWZm8uElm3zHX6v4A+zNU6mvZrwlM4pAsABApPB4uKRlO25ChzVbFprSIS6Jq+N3Mvrss+rnQvFTN+Y030Tz+Ksmt3PJHo0Ts9prlZW0ek+CweaJ1IlVzzN+9FlwX58FrfojLjewM6NSci7ZvpoRe59P+vi9J/23u3x5evOH1JCAAcECAJFAdnLxpF3Mm5Yc5k2buPrn39s5Py5dOlnJxq5XVjJ/9x3z3XczJzdQbplo3VrhddiBJ1Sl1LkKPoOeBPVqrT6x6hk9gNEQ7hsNMXWqer2RbyN5FC0eM0plNm/OHBvLD9ArHEOVTMT8888K219vnosaFjAgWADAavTkCwjUBRmfX36qiL9pMpbH0xw+n6RhmHfdVb3YbRuLuSI9UzqhBktLKz8J6sngqOdEqafMUNCTlwB5Fqp7dkOehYQE9f2plMExcFvLywzQHaKmPLfxBOV3Q+g8Bqv5NrNeOxAECwBEAj35AuQ6vRkUMzK4NK01L/77ch47lvn77/0l4hCxZvXO8P3nvssbmo5g7zs67nI9Hv0ZHPXeUekpMxT03Akig6Okc0sGx4kT1euXXBsbK7RwKB0zenVyjByDzNbVa4eBYAGASKH3rkN+J6tHq6GTDxGT3whPmCA89qj2yMKoX6N3UmbfeaFlIXpbFlq3Ft4+qWcbif1Igu2HzMygx4xvGfK/wYjiFgURI9fQWAIAhE5SkjGdHn1SkqauVSuiq68mSkiQpuXlEb34IlHPnkRt2xJNmkT0yy8h+tW7TqHq9ZQ3YQJR69ZEa9cSZWYS7dtHlJ1NlJ8vaI4eJTpzRvg/JYWoeXPh/7Nnhd/M1OTnC8vet0/wsnat4G3CBP3bM9j65OcTDRtGVFFBFB9PNGeO8LeiQpgu12VnEx04QJScTNSiheA5OZlo/37nlam0zRo00LeNDhwQytDYD8WZHalywj+DHjO+Zcj/6tlvQMLiwMVS0LIAajInTjDPncs8dGj1N/USMffubbfDMMFoiOgcDWHGNqoq87mny7hZM+b77mPOyzO+S2o6eAwBQA3jyBHmmTOZBw4UUksTMb/8sr/G6xVe23vggD0ew0LPOP5IaiKxPka9uKVME7dRly7ST/v3h76omgqCBQBqMIcOMU+fzvz77/7Tf/xROrH27cv8+uvs/4Irp4MMjsG9uKVMvWiUuXOnNPnCC8NfVE0EwQIAoBoPP1z9UUVMjPDOirffZi4stNuhBmhZQMtCQJlPPSVNfvXV0BdTk0GwAACoxt69zM8+y9y5c/WggUjo9zBkCPPHH9vtNAA9z68jqYnE+hjRualMk7aRN8/DbdtKAW9gKxrQB4IFAIAmO3YwT5rE3KZN9aBh7Fi73clQy6gnn66W8toKjbwToFXrk5mp/x0Gubn6yszIEMqUd97UWr/ATotaZerxqXebKXWsVChzS/MrffW1f//QdgVAsAAA0InXy7x5s5CjIT1dOPl+/bW/5tQp5jFjmD/7LMJDz9UuGiJTp/qP91dKeR2oiY8XphnVBAYVGRnGE/YEWx+1HALBdGJqZDWtmLZZHDIjT2WtVGZKSvDlB5YZzKeebVaVkEzPu04m0hTf7nrrxu+VywNBQbAAADCM18u8YQNzaan/9A8+kK6jdesy33QT81dfMZeUWGgm2IVVKZlSYLCglnApIyN4Uia5Ru5J/j4OJU2o6yOiN9NjYGKklBTlMouKhAyH8nVr0ULZt8dT/Z0MgYmWlMrUyvSod5sF7geNd5148zycEecRFk3l/EdqzxqRQMkKkJQJAGCYmBiiCy8kSkz0n750qfT/mTNEH34oJIVq2pTottuIvvlGyMljGsXFRIMGScl3Vq8mSk2trmOW/o+PJ/J4hPmKi5U1atP0aJTQoyHSvz7FxURjx0oJjyoqhO/y9RFp0ULYASKFhUQlJfr8qPkuKRHKEWnSRFiOEnFx0v+VlUSjRyv71LvsYL/J+HFbIh2oFLbfIFpFjWOP6poPhAeCBQCAJv/5jxAQjB9PdP750vSTJ4nee49o6FAh+eHUqSYtUJ7NT+3CmpRE9PjjROnpRGlp0gX29tv9M/Xdfbf/hS0+XphmREMkZQ/0eITlpacLyzeawVFtfQJ1wbJEJiURPfmkcDFXyqAocvSoFESI61haKmVFlK+fPDNj8+ZETz2lvmxx2ytlvJSXqXebBe5Pj8c/c6SsTM/t/6JkOkJERGMartC/H0B4WNzKYSl4DAFAZCktZV68WHgUUbeuf+v2M89U12u+pyIYJnaIc0QHR6veI2KkM6KeDo65uca2vd5Okybuz7KWWfzfD/7gE4fw+CEc0GcBAGA5RUXMCxYwjxrFnJRU/fH2L78YeMFVuLhp6KQVOHiYo2PKBNVAsAAAiChnz1af9uST/i0PbdoIwzV37LDIhJ7EQGZpnIiDEyg5pkzgBzo4AgAiSp061acVFvp3Bdizh+iZZ4g6dSLq3Jno2WeFaaaRmko0b57/tHnz/PsImKVxInp9W7F+ESizguLo6IyPnb8fohQECwAAS3j7baKCAqKZM4kGDhRGW4js2CG8Rjsri+iRR0xaYH4+0bhx/tPGjfPvJGeWxono9W3F+kWgzNWUTU2v7EGXDyym//0v9GJBaCBYAABYRnIy0V13Ea1aRXTwINH06UQXXeSv6dXL/3tpKdHhwzoXIA7XE3vei8MT160T/u7bR9S/v/B7fj7RgAGha+S98/UME4wEetY/O1tqwtmzR1tnZP30LtukMj859w6qpHhavjqJjv56VH+ZwBwsfiRiKeizAIA7yctjfukl5ksuYT5zxv+3zz+XXnA1c6bw+m1FZs1ibt1aO+WwmGkwJYW5efPqGQz1agJTMbdubTyDo9noWX95psc5c5QzM4ayfkaWbUKZpaXM9etVMBHzOXSaz2Z0cM5+cDHo4AgAcC2jR/v3aYuLYx46lPm995hPnKgSFRUJFwrxQqjU+U0tO2N6evAMjnINs3Iq5tat7cscqGf9matnelTLzGhk/fQu28QylyyR7P/lnC+dsx9cDjo4AgBcS7duRG3aSN8rK4n++18hKVSTJkTXXEP00edJdGbWR1JCovh4oo8+8u/8lpRENH++fy/L2FiiTz7xT8oUTEMklPtRwPI++MC+ZEBJScLytdafSNiQM2f6T5s5038DExlbP73LNrHM+fMl6Zin2jpnP9QkLA5cLAUtCwBEJ+ILrh59VHrBVeDnn/XfUrwTXbGCuX175hVfl6BlIQpaFpYsYY6NFSbVO6+SS1q2c85+cDloWQAAuJqYGKLu3YUU0vv3E+XkED3wgJBWWmT03xtLqZEzM+ngvhL66oJ/0cSHSmnXLqLHHiPi02cEcUqKcmrio0eFF15oaQI72wVLxRwJAlNDK3UoDEzhPGeOcmpoo+und9kmlMkDsun+v5WT1yvIRlQuolr7f3XOfqhBxDAz220iVAYPHkxERCtXrrTZCQAgEni9wjVi5Uqip58miikpFi4U+fn0QrePaOLxR/30PelHapR0lqhPH2HCxo10TfF8+r/MFcIY/nHjiPbtoyuTVpK3T1+fhoqLiZLqEHXtQrRtO1FxEVFSHfrXp+2o1/AmVb8n0ZYtQlCih6+/9n/aMWcO0cKFwefr1o3o+ef9p919N9GBAyQ8o4mLIyopJtr4g8+n6Pu24hl0feYWouXLidq0oeM//kY3Dfi9ms73vU9votpJRJWV9NY7cdSypbTMFSuIXnlFZkJl2fX7daCPfrvQ78VZUz9Mpe9WVfpvAAUuHVBOD09M8AUb/93XmobRf32/P0uP0+OZ86V3bFTtBxAaRq6h8VabAQAAs4iNFUY59u9fNUG8UKSm0vxmDxId99dvol5ExUS0WpySTVnnHSDa9w7RxRcLkzIzafmBgeRdHePTEJEwXw4R0QDf93vFttiq5R47JlyHQ2H3bn3zlpdXn7ZuHdHPPxMRiRffJD+fou9BDTYQrX7V9+y/tEUrWl7cqprO9321uIQ4On3af5kHDwb6VV520xWHicj/DZvbthEtX6EdKBARNWmSIPyTmkr8v9U0qc1xojImohgiYvo8cTQ99r//oxixfwQChYiBxxAAANfDTHTtDYnBhUREQ4f4f583j4SLURQyZox9GQ/DzOD4za5U+rGsK0n7JoY2lXWhb3Yhg6Md4DEEAMD1MBP16VZKm7fHUaWswTSOKqlLx0pa+X2ikEHy998p8aqhVOfATmnmzEw6+eV3wuueSdDQ8OFEeQckTXoG0ddf0zltW1BCgjS5vFzo3qCHevX8s1gWFwsJqIIRH09Ut67/tFOnyPcc34eC79oZzaj2mm98F22vl+jUTvX1820DIjr3XP+nBmVlREVFKiZly44hpnp0yq9l4exZ5RaSQBIThdTh6vuzgnp0qaSNW2v5bUsQGoauoRZ2tLQcjIYAADAzL59XqDhiQvwsn1eIt04a0VmxbAPo2p8gbJCUCQBQY/DmebhX4laOpQrFC0ssVXCvxK3sTU3TvnimpQmfYBdYpwUMav6Usijq0RlZP73LNlCm7v2Z57D94EIwdBIAUDPIz6ey7CHkKWtCXlLuQOelOMova0xl+QVEaWlST3oi4e/q1cJ0j0f4qGmUhgjaTeAQRS3f/frp0+ldPyPLNlCm7v2ZPcQ5+6EGgNEQAAB3UnWxqrV/H/2Ydj0deXsBUdOm/pqCAqLbbqPGBduoFpWFvizx4ideHLOz/S+OdqB1sRYRsyj26xc826KR9dO77BDK1NyfREQFBdT47uuFfAtO2A81BSuaNnJycjgrK0vxs23bNtX5xo4dW00/adIkVT0eQwBQQ5Fn/FNr5pZr5I8Y5Bn/9GjkyJvXnZLBUauZXykzopZvPeund9l2lwmCYvtjiO7du9PatWv9PqNGjaIWLVpQ586dNee94YYb/OabMGGCFRYBAG5GnvFP7c4yMCugUsY/PRo54t2y3ZkD9ax/oE5PxkM966d32XaXCUwlIkMny8vLqX///jR27Fi65557VHXjxo2jdu3a0eOPP66rXAydBKCGoyeDn1yjptejMbrcSKDXh6gzqjdj2XaXCVQxcg2NSAfHVatW0cmTJ2nkyJFBtYsXL6Y+ffrQ8OHD6eWXX6bi4uIIOAQAuBI9Fwq5Rk2vR2N0uZFArw95S4pZ5RrdBnaVCUwhIh0cFyxYQP369aOmSp1VZAwfPpxSUlKocePGtHv3bpo2bRrt37+f3njjjUjYBAAAAIAChoKFadOm0ezZszU1S5cupVatWvm+FxQU0Nq1a+m1114LWv7o0aN9/7dt25YaNWpEt9xyC3k8HkpLSzNiFQAAAAAmYShYGD9+PF177bWamtSAjikLFy6k888/nwYNGmTYXNeuXYmIKC8vD8ECAAAAYBOGgoUGDRpQgwYNdOuZmRYtWkQjRoygBHlCdZ3s2rWLiIgaNWpkeF4AAAAAmIOlHRxzcnLo999/p+uvv77ab3/88QcNGzaMtm/fTkREHo+H3nzzTdqxYwf9/vvvtHLlSnr00UepV69e1K5dOyttAgAAAEADSzs4LliwgLp37+7Xh0GkvLyc9u/f7xvtkJCQQBs2bKD333+fioqKqFmzZjRkyBD629/+ZqVFAAAAAAQBr6gGAAAAaiCOy7MAAAAAAPeCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgCYIFAAAAAGiCYAEAAAAAmiBYAAAAAIAmCBYAAAAAoAmCBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKAJggUAAAAAaIJgAQAAAACaIFgAAAAAgCYIFgAAAACgSbzdBsKhsLCQKisrafDgwXZbAQAAAFzF4cOHKS4uTpfW1S0LtWrVovh4V8c7AAAAgC3Ex8dTrVq1dGljmJkt9gMAAAAAF+PqlgUAAAAAWA+CBQAAAABogmABAAAAAJogWAAAAACAJggWAAAAAKCJ64OFt99+m8aMGUNdu3alCy64QFFz6NAhuvPOO6lr167Ut29feuGFF6iiokKz3JMnT9LDDz9MPXr0oAsuuIAee+wxOnv2rBWrYBobN26ktm3bKn62b9+uOt+4ceOq6Z988skIOjeHQYMGVVuPWbNmac5TWlpKkydPpj59+lD37t3p3nvvpaNHj0bIsTn8/vvv9Nhjj9GgQYOoS5cudOmll9L06dOprKxMcz437/cPP/yQBg0aRJ07d6ZRo0Zp1m8iomXLltGwYcOoc+fOdNVVV9F3330XIafm8c4779DIkSOpe/fu1LdvX/rb3/5G+/bt05xn0aJF1fZx586dI+TYXGbMmFFtXYYNG6Y5TzTsdyLlc1vbtm1p8uTJinor9rvrkxSUl5fTsGHDqFu3brRgwYJqv1dWVtJdd91FycnJNH/+fCosLKRHH32UEhIS6KGHHlIt9x//+AcdOXKE5s6dS+Xl5fTYY4/Rk08+SS+//LKVqxMW3bt3p7Vr1/pNe/3112nDhg1BK8oNN9xA9913n+97UlKSJR6t5r777qMbbrjB9/2cc87R1E+ZMoW+++47eu211+jcc8+lZ555hv7+97/T/PnzrbZqGvv27SNmpn/961+Unp5Oubm5NGnSJCouLqZHH31Uc1437velS5fS888/T5MnT6auXbvSf/7zH7rtttto+fLl1LBhw2r6zZs308MPP0wPPfQQDRw4kBYvXkz33HMPLVq0iLKysmxYg9D44Ycf6KabbqLOnTtTZWUlvfLKK3TbbbfRkiVLqE6dOqrz1a1bl5YvX+77HhMTEwm7ltCmTRuaO3eu77tWQqFo2e9ERAsWLKDKykrf9z179tCtt96qGSyZvt85Sli4cCH37Nmz2vTVq1dzu3bt+MiRI75pH330Effo0YNLS0sVy9q7dy9nZWXx9u3bfdO+++47btu2LRcUFJhv3iLKysr4wgsv5DfeeENTN3bsWH722Wcj5Mo6Bg4cyHPnztWtP3XqFHfs2JGXLVvmmybu+y1btphvMILMnj2bBw0apKlx636//vrrefLkyb7vlZWV3K9fP37nnXcU9ffffz/feeedftNGjRrFkyZNstSn1Rw7doyzsrL4hx9+UNWonRfdyPTp0/nqq6/WrY/W/c7M/Oyzz/Kll17KXq9X8Xcr9rvrH0MEY+vWrZSVlUXJycm+af369aMzZ87Q3r17FefZsmULnXfeeX534xdddBHFxsYGbe50EqtWraKTJ0/SyJEjg2oXL15Mffr0oeHDh9PLL79MxcXFEXBoPrNnz6Y+ffrQiBEjaM6cOZqPm3bs2EHl5eV00UUX+aa1atWKUlJSaOvWrRFwax2nT5+mevXqBdW5bb+XlZXRL7/84rfPYmNj6aKLLqItW7YozrN161bq27ev37R+/fpFxT4moqD7uaioiAYOHEgDBgygu+++m/bs2RMJe5aQl5dH/fr1o8GDB9PDDz9Mhw4dUtVG634vKyujr776ikaOHKnZWmD2fnf9Y4hgHD161C9QICLf9yNHjqjO06BBA79p8fHxVK9ePdV5nMiCBQuoX79+1LRpU03d8OHDKSUlhRo3bky7d++madOm0f79++mNN96IkFNzGDduHHXo0IHq1atHW7ZsoVdeeYWOHDlC//znPxX1R48epYSEBDrvvPP8pjds2NBV+zmQvLw8+uCDD4I+gnDjfj9x4gRVVlZWe9zQsGFD1ef3SueAhg0buq5vihyv10tTpkyhHj16aDapt2zZkqZMmUJt27al06dP03vvvUdjxoyhJUuWBD0vOI0uXbrQ888/Ty1btqQjR47Qm2++STfddBMtXryY6tatW00fjfudiOjbb7+l06dP07XXXquqsWK/OzJYmDZtGs2ePVtTs3TpUmrVqlWEHNlLKNujoKCA1q5dS6+99lrQ8kePHu37v23bttSoUSO65ZZbyOPxUFpaWsi+zcDIut96662+ae3ataOEhAR66qmn6OGHH6bExESrrZpOKPv9jz/+oNtvv52GDRvm13dDCSfvd6DN5MmTac+ePfTRRx9p6rp3707du3f3+37FFVfQ/Pnz6YEHHrDYpbkMGDDA93+7du2oa9euNHDgQFq2bBmNGjXKRmeRZeHChdS/f39q0qSJqsaK/e7IYGH8+PGaURMRUWpqqq6ykpOTqz06ECPLRo0aqc5z/Phxv2kVFRX0559/qs5jJaFsj4ULF9L5559PgwYNMry8rl27EpFwh2r3RSOcutC1a1eqqKig33//nTIzM6v9npycTOXl5XTq1Cm/1oVjx47Zsp8DMbruf/zxB/31r3+l7t270zPPPGN4eU7a72rUr1+f4uLi6NixY37Tjx07Vu0uUiQ5Obna3aSW3un861//otWrV9MHH3xg+C4xISGB2rdvTx6PxyJ3keO8886jjIwM1XWJtv1ORHTw4EFav349zZgxw9B8Zux3RwYLDRo0qPYYIFS6detGM2fOpGPHjvmaLtevX09169al1q1bK87TvXt3OnXqFO3YsYM6depEREQ5OTnk9XqpS5cupvgygtHtwcy0aNEiGjFiBCUkJBhe3q5du4hIPZiKJOHUhV27dlFsbKxiD3kiok6dOlFCQgJt2LCBhg4dSkTCyIJDhw5Rt27dQrVsGkbWXQwUOnbsSM8//zzFxhrvjuSk/a5GYmIidezYkTZs2ECXXnopEQlN8hs2bKCxY8cqztOtWzfKycmhW265xTdt/fr1jtjHRmBmeuaZZ2jFihU0b9483TdMciorKyk3N9fvLt2tnD17lvLz81Xra7TsdzmLFi2ihg0bUnZ2tqH5TNnvpnaXtIGDBw/yzp07ecaMGdytWzfeuXMn79y5k8+cOcPMzBUVFTx8+HAeP34879q1i9esWcMXXnghv/zyy74ytm3bxkOHDvUb6XDbbbfxiBEjeNu2bfzTTz/xkCFD+KGHHor4+oXC+vXrOSsri/fu3Vvtt4KCAh46dChv27aNmZnz8vL4jTfe4J9//pnz8/P522+/5cGDB/NNN90UadthsXnzZp47dy7v2rWLPR4Pf/nll3zhhRfyhAkTfJrAdWdmfvLJJzk7O5s3bNjAP//8M48ePZpHjx5txyqETEFBAV922WV88803c0FBARcWFvo+ck207PclS5Zwp06deNGiRbx3716eNGkSX3DBBb4RT4888ghPmzbNp9+0aRN36NCB3333Xd67dy9Pnz6dO3bsyLt377ZrFULiqaee4p49e/LGjRv99nFxcbFPE7juM2bM4O+//549Hg/v2LGDH3zwQe7cuTPv2bPHjlUIi6lTp/LGjRs5Pz+fN23axLfccgv36dOHjx07xszRu99FKisrOTs7m1966aVqv0VivzuyZcEI06dPp88//9z3fcSIEURE9P7771OfPn0oLi6OZs6cSU8//TSNHj2akpKS6Nprr/UbW15cXEz79++n8vJy37Rp06bRM888QzfffDPFxsbSkCFD6IknnojYeoXDggULqHv37op9OsrLy2n//v2+Xu/infX7779PRUVF1KxZMxoyZAj97W9/i7TtsEhMTKSlS5fSG2+8QWVlZdSiRQu65ZZb/PoxBK47EdFjjz1GsbGxdN9991FZWRn169ePnnrqKTtWIWTWrVtHeXl5lJeXR/379/f7bffu3UQUXfv9iiuuoOPHj9P06dPpyJEj1L59e5ozZ46vefnw4cN+LSs9evSgadOm0WuvvUavvPIKZWRk0Jtvvum6sfYff/wxEQkdeeU8//zzdN111xFR9XU/deoUTZo0iY4cOUL16tWjjh070vz581VbVZ1MQUEBPfTQQ3Ty5Elq0KAB9ezZkz799FNf61u07neR9evX06FDhxRHt0Viv8cwM4c8NwAAAACinqjPswAAAACA8ECwAAAAAABNECwAAAAAQBMECwAAAADQBMECAAAAADRBsAAAAAAATRAsAAAAAEATBAsAAAAA0ATBAgAAAAA0QbAAAAAAAE0QLAAAAABAk/8Hz4A7NEMV2EIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use('seaborn-white')\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.axis('equal')\n",
    "for idx in range(len(views)):\n",
    "    color = 'r'\n",
    "    cam_center = cam_center_list[idx]\n",
    "    plt.scatter(cam_center[0], cam_center[dim], s=100, c=color, marker='x')\n",
    "    # plt.text(cam_center[0], cam_center[1], f'{idx}', fontsize=12)\n",
    "\n",
    "cam_center_array = np.array(cam_center_list)\n",
    "hull = ConvexHull(np.array(cam_center_array)[:, [0, dim]], 'Qg')\n",
    "hull_list=hull.vertices.tolist()\n",
    "\n",
    "simplified_hull = []\n",
    "dist = 3\n",
    "for idx in range(len(hull_list)):\n",
    "    dist_1 = np.linalg.norm(cam_center_array[hull_list[idx], [0, dim]] - cam_center_array[hull_list[idx-1], [0, dim]])\n",
    "    if dist_1 > dist:\n",
    "        simplified_hull.append(hull_list[idx])\n",
    "        \n",
    "simplified_hull = simplified_hull[2:] + simplified_hull[:2]\n",
    "simplified_hull.append(simplified_hull[0])\n",
    "\n",
    "plt.plot(cam_center_array[simplified_hull,0], cam_center_array[simplified_hull, dim], 'b--^',lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.2 Plot Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "with torch.no_grad():\n",
    "    voxel_size = torch.tensor([0.25, 0.25])\n",
    "    xy_range = torch.tensor([-9, -7, 7.5, 6])\n",
    "    xyz = gaussians.get_xyz.cpu()\n",
    "    voxel_index = torch.div(torch.tensor(xyz[:, :2]).float() - xy_range[None, :2], voxel_size[None, :], rounding_mode='floor')\n",
    "    voxel_coords = voxel_index * voxel_size[None, :] + xy_range[None, :2] + voxel_size[None, :] / 2\n",
    "    voxel_dim = torch.tensor([int((xy_range[2] - xy_range[0]) / voxel_size[0]), int((xy_range[3] - xy_range[1]) / voxel_size[1])])\n",
    "    print(f\"Dimension of Voxels: {voxel_dim}\")\n",
    "\n",
    "    new_coors, unq_inv, unq_cnt = torch.unique(voxel_coords, return_inverse=True, return_counts=True, dim=0)\n",
    "    feat_mean = torch_scatter.scatter(xyz[:, -1], unq_inv, dim=0, reduce='mean')\n",
    "    feat_std = torch_scatter.scatter_std(xyz[:, -1], unq_inv, dim=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAGuCAYAAADf4Q+iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8MklEQVR4nO3deXwU9eH/8XcSCOESNUFBQBAEhCRcoiYYjBA0eKBYA3KJVUTxqEJB0CpV+AJBESkgVhSECCilglRRURuFarkCQi385KwHCEgCyBUkkMzvj3nMHslusrtkZ3fh9Xw89rHZ2c/MfOYzM59578zsJsowDEMAAAA2iA51BQAAwPmD4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAUSwli1bavr06SGZd9euXfX000+HZN7evPnmm+revbtKSkpCXZWgOXz4sNq1a6eVK1eGuipAQAgegBdDhgxR27Ztdfz4ca9lhg8frqSkJB0+fNjGmvlnz549atmypWbPnh3qqgTV8ePHNWvWLA0ePFjR0dFl3nvppZfUtWtXJSUlqXPnznriiSd08uRJR5nVq1frmWeeUWZmptq2bauMjAw9++yzOnDgwFnV669//auGDBmiTp06lRsUP/vsMw0dOlQZGRlq27atMjMzNXHiRB09etSt3EUXXaSsrCxNnTr1rOoFhEqVUFcACFd33HGHvvzyS/3zn/9Uz549y7x/8uRJffHFF0pLS9NFF11kfwXh5r333tOZM2d0++23uw0/duyYBgwYoP379+uee+7R5ZdfrkOHDmnDhg0qKipS9erVJUmTJk3SkSNH1L17dzVp0kS7d+/W/PnztWLFCi1dulR169YNqF5/+ctfVLduXbVq1Upff/2113KjR4/WJZdcojvuuEOXXXaZtm3bpvnz52vlypV6//33FRcX5yjbt29fzZs3T6tXr1ZqampA9QJCheABeNG1a1fVrFlTH374ocfgkZubq8LCQt1xxx32Vw5lLFmyRF27dlW1atXchk+ePFl79+7VkiVL1KhRI6/jP/PMM7r66qvdzpZ07txZAwYM0Pz58zVs2LCA6pWbm6uGDRvq0KFD5YaEadOm6brrrnMblpSUpFGjRunDDz9Ur169HMObNWumFi1a6P333yd4IOJwqQXwIi4uTjfffLPWrFmjgwcPlnl/2bJlqlmzprp27SpJ2r17t5544glde+21atu2rXr37q0VK1Y4yu/atUtt2rTRyJEj3aazfv16tWrVSpMmTXIMO3r0qMaPH6/09HQlJSXppptu0htvvFFp9y4sWbJELVu21IYNG5Sdna2UlBS1a9dOjz32mA4dOuRW1jAMvfbaa7rhhhvUtm1b3XvvvdqxY4fH6VZUb8MwdO+99yolJcWtTYuKitSjRw9169ZNhYWFjuG7du3S3r17K1ye3bt3a9u2berUqVOZ+ixZskS9e/dWo0aNVFRUpKKiIo/TuOaaa8pcornmmmt04YUX6n//+1+FdfCmYcOGPpUrHTokqVu3bpLMdiitU6dO+vLLL8U/GEekIXgA5ejRo4fOnDmjTz75xG34r7/+qq+//lo33XST4uLiVFBQoD59+ujrr79W3759NWzYMJ06dUqPPPKIPv/8c0nmp9Qnn3xS//jHP5SbmytJKiws1DPPPKOmTZvqySeflGRewhkwYIA++OAD9ezZU88995w6dOigV155RdnZ2ZW6fOPGjdPWrVv1+OOPq2/fvvryyy81duxYtzJTp07V1KlTddVVV2nkyJFq1KiRHnjgAbeA4Gu9o6KiNGHCBJ06dUrPP/+8Y9zp06drx44dys7OVo0aNRzDb731Vo0aNarC5di4caMkqXXr1m7DN2zYoFOnTqlx48Z64okn1K5dO7Vp00Z9+vTRd999V+F0T5w4oRMnToTsUlpBQYEkeZx/YmKijh496jUEAuGKSy1AOVJSUlS3bl0tW7ZMAwYMcAxfvny5Tp8+rR49ekiS3njjDRUUFGjBggXq2LGjJKlXr1664447lJ2drYyMDEVHR+v+++9Xbm6u/vznP6tDhw6aPn269u7dq4ULFyo2NlaSNGfOHO3evVvvv/++mjRpIknq06ePLrnkEs2ePVsPPPCA6tevXynLd+GFF+qtt95SVFSUJKmkpETz5s3TsWPHVLt2bR06dEizZs3SjTfeqNdff91RbsqUKXr99dfdpuVrvRs1aqSnn35af/7zn/XBBx+ocePGmj17tgYOHKhrrrkmoOWwzkiUPrvw448/SjIvt1x++eV68cUXdezYMc2YMUP33Xefli1bpksuucTrdHNycnT69GndcsstAdXrbL355puKiYlRZmZmmfesy0Y7d+5UixYt7K4aEDDOeADliImJ0W233aaNGzdqz549juHLli1TQkKC4/r6ypUr1aZNG0fokKSaNWvqnnvu0c8//6ydO3dKkqKjozVx4kQVFhZq8ODBeuedd/TQQw8pOTnZMd7y5ct19dVX64ILLtChQ4ccj06dOqm4uFh5eXmVtny9e/d2hAlJ6tixo4qLi/Xzzz9LklatWqXTp09rwIABbuXuu+++MtPyp9733HOP0tLSNG7cOMdZlD/+8Y9lprlt2zbNmzevwuX49ddfVaVKFdWsWdNt+IkTJySZZ1rmzp2rHj16qF+/fpoxY4aOHDmiBQsWeJ1mXl6eZsyYoVtuuSUk91F8+OGHeu+993T//fc7gpyrCy64QJLC+htVgCec8QAq0KNHD82dO1fLli3TkCFDtH//fq1fv1733nuvYmJiJEl79+5V27Zty4zbtGlTx/vWp9LLL79cjz/+uF566SW1aNFCjz76qNs4P/74o7Zt2+b1YFf6Hoyzcdlll7m9tg5m1lc4rfsrSh/4Lr74YtWpU8dtmL/1njBhgrp166YjR45o4cKFbt/aqCzWNLt06eIWStq1a6eGDRs6LtGUtmvXLj3++ONq3ry5xo0bV+n1qsj69ev17LPPKi0trcKbWl0DIRAJCB5ABZKSktS0aVN99NFHGjJkiJYtWybDMByXWQLx73//W5J04MAB/frrr25f1SwpKdH111+vBx980OO4nj79Bqr0zZSWQG5Y9Lfea9euddzouX37drVv397veVouvPBCnTlzRsePH1etWrUcw63LKAkJCWXGiY+PL/MbGZK0b98+DRo0SLVq1dIbb7zhNj07bN26VY888oiaN2+uadOmqUoVz930kSNHJHm+/wMIZwQPwAc9evTQ1KlTtXXrVi1btkxNmjRRmzZtHO9fdtll+v7778uMZ9174Hpm4d1339W///1vDRs2TDNnztSf//xn/fWvf3W8f/nll6uwsLDMNzRCwar3Dz/84PZV1EOHDjkOfBZ/6n3gwAGNGzdOaWlpqlq1ql588UWlpaWpQYMGAdXTOrO0Z88eXXXVVY7hiYmJkqRffvnFYx2s8SyHDx/WAw88oKKiIr3zzjvl3v8RDD/99JMefPBBXXzxxXrzzTfLXDpyZV36a9asmV3VAyoF93gAPrDObkybNk3fffddmbMd6enp+vbbb91O3RcWFmrRokVq0KCBrrzySknm1z5feuklZWZmasiQIRo1apS++OILLV261DHeLbfcoo0bN+qrr74qU4+jR4/qzJkzQVhCzzp16qSqVatq/vz5bmdBcnJyypT1p96jR49WSUmJxo8fr7Fjx6pKlSp69tlny5xp8fXrtNbZks2bN7sNb9q0qa666irl5ua6Xer5+uuvtW/fPreQVFhYqIceeki//PKL3njjjUo9s+SL/Px8PfDAA4qKitLs2bN18cUXl1t+y5Ytql27tpo3b25TDYHKwRkPwAeNGjVS+/btHV+DLR08HnroIX300UcaPHiw7r33XtWpU0dLly7Vnj17NH36dEVHR8swDP3pT39SXFycXnjhBUnmtz4+++wzjR8/Xqmpqbr00ks1aNAgffHFFxoyZIjuuusuJSYm6uTJk9q+fbs+/fRT5ebmVnhQqiwXX3yxHnjgAc2cOVMPP/yw0tPT9f/+3//Tv/71rzKn+H2t9+LFi7VixQpNnDhR9erVkyQ999xzeuqpp/TOO++of//+jmneeuutuvbaayu8wbRRo0Zq0aKFVq9eraysLLf3nnnmGT3wwAPq16+f+vTpo2PHjmnOnDlq0qSJ+vbt6yg3YsQIffvtt7r77ru1a9cut9/OqFmzpuM3NSTz67+vvvqq3n77bY+/v+Fq6dKl2rt3r3777TdJ5k2rr732miTpzjvvdJzlefDBB7V79249+OCD2rBhgzZs2OCYRkJCgq6//nq36a5atUpdunThHg9EHIIH4KMePXpo48aNatOmjRo3buz2XkJCghYuXKhJkyZp/vz5OnXqlFq2bKnXX39dN954oyRp3rx5WrdunaZPn+4WHMaPH6/bb79do0eP1htvvKHq1atr3rx5mjlzppYvX66lS5eqVq1aatKkif7whz+odu3adi62hg4dqtjYWC1cuFBr165VmzZt9NZbb+nhhx92K+dLvffv36/s7Gx16dJFd911l2PcO+64Q5999plefvll3XDDDeX+wqg3d999t6ZOnarffvvN7UbVlJQUzZo1S1OnTtUrr7yi6tWrq1u3bnrqqafcLmVs3bpVkrR48WItXrzYbdoNGjRwCx6FhYWKioryeO9IaYsXL9a6descr9euXau1a9dKkq6++mpH8LDmP2vWrDLTuPbaa92Cx65du7R9+3b96U9/qnD+QLiJMvjZOwDngGPHjqlbt24aMWKE28+LB0NWVpYuu+wyTZs2Lajz8Wb8+PFav369lixZwhkPRBzu8QBwTqhdu7YGDRqk2bNnV9pPy3ty/Phxbd261fFLs3Y7fPiw3nvvPQ0dOpTQgYjEGQ8AAGAbzngAAADbEDwAAIBtCB4AAMA2BA8AAGCbsPsdj44dO6qoqMjtf1cAAIDwlp+fr9jYWK1fv77ccmEXPE6dOqXi4uJQVwMAAPjhzJkzPv2DybALHtY/ZbJ+mhoAAIS/jIwMn8pxjwcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbONX8HjnnXfUo0cPdejQQR06dNA999yjlStXljvOJ598ou7duys5OVk9evSosDwAADh3+RU86tWrpxEjRmjJkiVavHixUlJS9Nhjj2nHjh0ey3/zzTcaPny4srKytHTpUmVkZOixxx7T9u3bK6XyAAAgsvgVPLp27ar09HQ1adJEV1xxhYYNG6YaNWpo06ZNHsu//fbb6ty5sx588EE1a9ZMQ4cOVevWrTV//vzKqDsAAL5Zv17q3998RkgFfI9HcXGxPvroIxUWFqp9+/Yey2zatEmpqaluw9LS0rwGFQAAgmLKFGnRIvMZIVXF3xG2bdumPn366NSpU6pRo4ZmzJihK6+80mPZgoICJSQkuA2Lj49XQUFBYLUFACAQw4a5PyNk/A4eV1xxhZYuXapjx47p008/1ahRozR//nyv4QMAgJDr2FFasCDUtYACuNQSGxurxo0bKykpScOHD9dVV12lt99+22PZhISEMmc3Dh48WOYsCAAAOD+c9e94lJSUqKioyON77dq105o1a9yGrVq1Su3atTvb2QIAgAjkV/CYPHmy8vLytGfPHm3btk2TJ0/WunXr1KNHD0nSyJEjNXnyZEf5gQMH6quvvtJbb72lXbt2afr06dq8ebMGDBhQuUsBAAAigl/3eBw8eFCjRo3SgQMHVLt2bbVs2VKzZ8/W9ddfL0nat2+foqOdWaZDhw56+eWX9Ze//EWvvPKKmjRpohkzZqhFixaVuxQAACAiRBmGYYS6Eq4yMjIkSbm5uSGuCQAA8JWvx2/+VwsAALANwQMAANiG4AEAOD/xM+oh4fcPiAEAcE6wfkZd4sfFbETwAACcn/gZ9ZAgeAAAzk/8jHpIcI8HAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANimij+FZ86cqc8++0z/+9//FBcXp/bt22vEiBFq2rSp13GWLFmiZ555xm1YbGys/vvf/wZWYwAAELH8Ch7r1q1T//79lZycrOLiYr3yyisaNGiQPvroI9WoUcPreLVq1dLy5csdr6OiogKvMQAAiFh+BY/Zs2e7vZ44caJSU1O1ZcsWXXPNNV7Hi4qKUt26dQOrIQAAOGf4FTxKO3bsmCSpTp065ZYrLCxUly5dVFJSotatW+uPf/yjmjdvfjazBgAAESjgm0tLSko0YcIEdejQQS1atPBa7oorrtCECRP02muvadKkSTIMQ3369NH+/fsDnTUAAIhQAZ/xGDNmjHbs2KF33nmn3HLt27dX+/bt3V7feuutWrhwoYYOHRro7AEAQAQKKHiMHTtWK1as0Pz581WvXj2/xq1atapatWqln376KZBZAwCACObXpRbDMDR27Fh9/vnnysnJUaNGjfyeYXFxsbZv387NpgAAnIf8OuMxZswYLVu2TK+99ppq1qyp/Px8SVLt2rUVFxcnSRo5cqQuvfRSDR8+XJL06quvql27dmrcuLGOHj2q2bNna+/everVq1clLwoAAAh3fgWPd999V5J07733ug3Pzs7W7373O0nSvn37FB3tPJFy9OhRjR49Wvn5+apTp44SExO1cOFCXXnllWdbdwAAEGGiDMMwQl0JVxkZGZKk3NzcENcEAAD4ytfjN/+rBQAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALCNX8Fj5syZuvvuu9W+fXulpqbq0Ucf1f/+978Kx/vkk0/UvXt3JScnq0ePHlq5cmXAFQYAAJHLr+Cxbt069e/fX4sWLdKcOXN05swZDRo0SIWFhV7H+eabbzR8+HBlZWVp6dKlysjI0GOPPabt27efdeUBAEBkiTIMwwh05EOHDik1NVXz58/XNddc47HM0KFDdfLkSc2cOdMxrHfv3rrqqqs0duzYMuUzMjIkSbm5uYFWCwAA2MzX4/dZ3eNx7NgxSVKdOnW8ltm0aZNSU1PdhqWlpWnTpk1nM2sAABCBAg4eJSUlmjBhgjp06KAWLVp4LVdQUKCEhAS3YfHx8SooKAh01gAAIEJVCXTEMWPGaMeOHXrnnXcqsz4AAOAcFtAZj7Fjx2rFihXKyclRvXr1yi2bkJBQ5uzGwYMHy5wFgQ3Wr5f69zefAQAIAb+Ch2EYGjt2rD7//HPl5OSoUaNGFY7Trl07rVmzxm3YqlWr1K5dO78qikowZYq0aJH5DABACPgVPMaMGaMPPvhAkydPVs2aNZWfn6/8/Hz99ttvjjIjR47U5MmTHa8HDhyor776Sm+99ZZ27dql6dOna/PmzRowYEDlLQV8M2yY1Lu3+QwAQAj4dY/Hu+++K0m699573YZnZ2frd7/7nSRp3759io525pkOHTro5Zdf1l/+8he98soratKkiWbMmFHuDakIko4dpQULQl0LAMB57Kx+xyMY+B0PAAAijy2/4wEAAOAPggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG7+DR15enoYMGaK0tDS1bNlS//znP8stv3btWrVs2bLMIz8/P+BKAwCAyFTF3xEKCwvVsmVL3X333Xr88cd9Hm/58uWqVauW43V8fLy/swYAABHO7+CRnp6u9PR0v2cUHx+vCy64wO/xAADAucPv4BGonj17qqioSM2bN9fjjz+uq6++2q5ZAwCAMBH04FG3bl2NGTNGSUlJKioq0t///ncNHDhQixYtUmJiYrBnDwAAwkjQg0fTpk3VtGlTx+sOHTpo9+7dmjt3riZNmhTs2QMAgDASkq/TJicn66effgrFrAEAQAiFJHhs3bpVdevWDcWsAQBACPl9qeXEiRNuZyv27Nmj7777TnXq1NFll12myZMn65dfftFLL70kSZo7d64aNmyo5s2b69SpU/r73/+uNWvW6K233qq8pQAAABHB7+CxefNmDRw40PE6OztbknTXXXdp4sSJys/P1759+xzvnz59Wi+++KJ++eUXVa9eXS1atNCcOXOUkpJSCdUHAACRJMowDCPUlXCVkZEhScrNzQ1xTQAAgK98PX7zv1oAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2/gdPPLy8jRkyBClpaWpZcuW+uc//1nhOGvXrtVdd92lpKQk3XTTTVqyZElAlQUAAJHN7+BRWFioli1b6vnnn/ep/O7du/Xwww/ruuuu0z/+8Q/dd999eu655/TVV1/5XVkAABDZqvg7Qnp6utLT030uv3DhQjVs2FBPP/20JKlZs2basGGD5s6dq86dO/s7ewAAEMGCfo/Hpk2blJqa6jYsLS1NmzZtCvasAQBAmAl68CgoKFBCQoLbsISEBB0/fly//fZbsGcPAADCCN9qAQAAtgl68EhISFBBQYHbsIKCAtWqVUtxcXHBnj0AAAgjQQ8e7dq105o1a9yGrVq1Su3atQv2rAEAQJjxO3icOHFC3333nb777jtJ0p49e/Tdd99p7969kqTJkydr5MiRjvJ9+vTR7t279dJLL2nXrl1asGCBPvnkE/3+97+vnCUAAAARw++v027evFkDBw50vM7OzpYk3XXXXZo4caLy8/O1b98+x/uNGjXSzJkzlZ2drbffflv16tXTuHHj+CotAADnoSjDMIxQV8JVRkaGJCk3NzfENQEAAL7y9fjNt1oAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADYhuABAABsQ/AAAAC2IXgAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2wQUPBYsWKCuXbsqOTlZvXr10rfffuu17JIlS9SyZUu3R3JycsAVBgAAkauKvyN8/PHHys7O1pgxY9S2bVvl5ORo0KBBWr58ueLj4z2OU6tWLS1fvtzxOioqKvAaAwCAiOX3GY85c+aod+/euvvuu3XllVdqzJgxiouL0+LFi72OExUVpbp16zoeCQkJZ1VpAAAQmfwKHkVFRdqyZYs6derknEB0tDp16qSNGzd6Ha+wsFBdunRRenq6HnnkEe3YsSPwGgMAgIjlV/A4fPiwiouLy1xSiY+PV0FBgcdxrrjiCk2YMEGvvfaaJk2aJMMw1KdPH+3fvz/wWgMAgIjk9z0e/mrfvr3at2/v9vrWW2/VwoULNXTo0GDPHgAAhBG/znhcdNFFiomJ0cGDB92GHzx40Of7NqpWrapWrVrpp59+8mfWAADgHOBX8IiNjVViYqJWr17tGFZSUqLVq1e7ndUoT3FxsbZv3666dev6V1MAABDx/L7Ucv/992vUqFFKSkpSmzZtlJOTo5MnT+p3v/udJGnkyJG69NJLNXz4cEnSq6++qnbt2qlx48Y6evSoZs+erb1796pXr16VuyQAACDs+R08br31Vh06dEjTpk1Tfn6+WrVqpVmzZjkutezbt0/R0c4TKUePHtXo0aOVn5+vOnXqKDExUQsXLtSVV15ZeUsBAAAiQpRhGEaoK+EqIyNDkpSbmxvimgAAAF/5evzmf7UAAADbEDwAAIBtCB4AAMA2BA8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheLhav17q3998BgAAlc7vn0w/p02ZIi1aZP69YEFo6wIAwDmI4OFq2DD3ZwAAUKkIHq46duRMBwAAQcQ9HgAAwDbn1RmPEye8vxcTI8XF+VY2OlqqXj2wsoWFkrf/BxwVJdWoEVjZkyelkhLv9ahZM7Cyv/0mFRdXTtkaNcx6S9KpU9KZM5VTtnp1s50lqahIOn26csrGxZnbhb9lT582y3tTrZpUpYr/Zc+cMdvCm9hYqWpV/8sWF5vrzpuqVc3y/pYtKTG3tcooW6WK2RaSuU8UFlZOWX/2e/oIz2XpI/wvGw59REgZYaZr165G165dgzJtcxf1/Lj1VveyNWp4L5ue7l42IcF72Y4d3cs2buy9bOvW7mVbt/ZetnFj97IdO3ovm5DgXjY93XvZGjXcy956a/nt5iorq/yyx487y953X/llDxxwln300fLLfv+9s+yIEeWX3bzZWfb558svu26ds+xLL5Vf9ssvnWVffbX8ssuWOcvOmVN+2UWLnGUXLSq/7Jw5zrLLlpVf9tVXnWW//LL8si+95Cy7bl35ZZ9/3ll28+byy44Y4Sz7/ffll330UWfZAwfKL3vffc6yx4+XXzYry3BTXln6CPNBH+F8RHIfEQy+Hr+51AIAAGwTZRiGEepKuMrIyJAk5ebmVvq0OY3qf1lOo/pfNhxOo3KpxbeyXGpxoo/wv2wk9xHB4Ovx+7wKHgAAIDh8PX5zqQUAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwCLb166X+/c1nAADOc1VCXYFz3pQp0qJF5t8LFoS2LgAAhBjBI9iGDXN/BgDgPMallmDr2NE809GxY6hrAiCccBkW5ymCBwCEgnUZdsqUUNckuAhYoRWG7U/wQGDCcGNGkLCug2PYMKl373P/MuyUKdLf/iZlZVX+NhSKbTPS9ocwDLgEDwQmDDdm20Rax3O2zuV1vX69lJoqXXaZ9Pbb9s53yhQzdJzrl2GHDZMaNpT27Kn8bSgU26Y/83z7balJE/u3Ldf+KRwDrhFmunbtanTt2jXU1QiOvDzD6NfPfI5059Ky+KtfP8OoUsV8PhdUtC7P5XXdr59hSOajcePKmaYv7VV6G/JlnEheD8GqeyjaxJ95Nm7s+7bl77Lk5JjTzclxHx7C/snX4zfBw07n2gHLbuHS8dpVD7vmY/d2GS7r0apLUpJhxMUZxgsvVM40rfZMSSl7YMjLM4zMTPO91FTzb2tYVJTzdb9+5niuzykpZpl69cKj7VAxb+HAE1/2Q9d9p359M9TUr+8+PCfHHJaS4nztax3OEsEjHIVTh2sY4VefipwPwc11nbgubzDXld3bQajWo7flDKQ+5bWZ1fHHxjo/7VrlrYARFWUOr1LFHFavnjN4WPWx3reerQONZJYLlUjrN0LJn7by90xZSoq5LcTGmiHWGt6vn3Mb69fPv7MuZ4ngEa4qc6c922lF2oH8fOjwvIWNSFtX5QnVaXdvbehrfVw/OZa3PqyOXzKMGjXMMymNGxtGTIwZGDIzzQNFSorzdUyMM6B4O+NhBRrX4BGKfeJc2hb9EUhbV3ZbudYhL8/cviQzuFrbUk6O82/OePgmJMHDOtVprahAp+HLRlnRhujPxn22G/X5cCD3Rzi0h7c6hEPdwl1F902U14a+tK/rJ8fy1lNmpmEkJ5sHAys0uAaL0nX1p/8pPV9/+4DK2I7O122xss+MBaL09F54wRluzzZYVwKChzeeVkLpU1OB8LTSPc0r0E9lvi5LRfOwu9OIpE4qnD7JRVK7hYuzOSi7hgPrDEXpT4i+fHIsfZnE22Uyb3X1dF+IP8vsy3IGYxsPp0uBgdalovGCGSIC/eDqy+VYG/u1oAaP+fPnG126dDGSkpKMrKws4z//+U+55T/++GMjMzPTSEpKMm6//XZjxYoVXssGPXh4CwjePnEEchrWGs86vWp3Qna9Ua30tP2p09nUxfWadmVs9L6ErNKnqCvaof35RFyZXLcVfzqLUF2iCPb8y5t+6fVr7afewkHp6aWkmJcnKjqQu+4b1ulr12vivoZ8a91an0Bd6136hlJP03G9YbD0tnw2bWe9du3nrPd9acuK5lXR9lrR3+Xts57OGJUXAsvr/0rX29fLZ/60ha/DrPllZrr3y6W3IddLbdZ2ZF2ms9bdCy94Xrf+hJpKELTg8dFHHxmJiYnGe++9Z+zYscN47rnnjI4dOxoFBQUey2/YsMFo1aqV8eabbxo7d+40pkyZYiQmJhrbtm07q4r7zdrprB0/J8c8HVq1qmE0a2b+XbOmYVSrZhg9ezpXZv365kacmmo+4uPNstaKt1ZmgwZmhxEXZ067Xj3ztVW+Zk3DiI42p1W1qvsd9Hl5Ze+sdx02eHDZTqxZM/O9nj2ddzBbG2izZua869Rxbqypqc4b2GJizPl468hfeMGso3WdukEDsy4xMeYy1Knj3rHm5JjLWaeOWY/Bg53jRkebNz9Zy+V6x/XgweZ7Vj1d15PrtUqrHawb8lx3xMaNzeGxseZOLJl1t+Zt1SM+3myX2FjzOSXFnK9kGLVqme+npprTTEkxX1ttWro+rgePwYPNdklPN6cXHW2+Tkpy7xRc28q6QTAmxv1egMGDndtJerr7nemG4byZzBrmqcMu75S9p/dzcpwHWmtdWO1QevrWAaB+ffdp5OWZbRodbS6D6/xSU81puk7XUx3q1zfXsWROyyrjev+DtfzWfmYdoL11uIbh+fLI4MHO9vYUPF3ft7aT2Fj37c+6nFKtmvtwKzikpDjbLz7eWW/r0ayZuT82a2a+36CBud3GxpoPqz+pWdO5/bquE2sbTU01+wBPN50mJTn7AWvdW+s5Odm5/Vv7TZ06ZetUrZr5OjXVnF7Nmu77VGqqYdSta76uXdt9m46ONp+jopzziI01y5Te/q1yjRs7pycZxiWXOPdnq0yDBs5pl37ExZl1di1fvbqzPa1HzZrm8ljzr1rVeVBPTjaX7aKLnO1Tevtt1swcbo1v9YnW+pLM6aekuNe/Vi1z+lYbuLZTMB7R0YaRkODs+6xhSUlBCSJBCx5ZWVnGmDFjHK+Li4uNtLQ0Y+bMmR7LP/nkk8ZDDz3kNqxXr17G6NGjPZYPWvCwdjprQ6pVq3JXbnmvvT1cNz7XYa4bqj8Pa+f2pS6eygbyqFrV83J4enTo4N7RlH40aODsDMPx4dqpSOUvi6dHoNucFRZd29z1IBwVZR7wXA9wsbHm/GJinKGsdH1Ld/zlPfzpHK2gYYUDb4+YGLOO5dXBNTz6Wy/X8WrX9l62enXP+5y3+QbzQBGsR2xs4P3K+fTwd5+O5EdqaqUfZn09fvv1y6VFRUXasmWLOnXq5BgWHR2tTp06aePGjR7H2bRpk1JTU92GpaWladOmTf7/2lllOXNGOn688qZXUlL+a2+Kiz1P6/TpwOpx5ozn6UlSVFTFZQNx+rTn5fDkm2/MTd6bn3+Wjhzx/n7pZbDbiRPur8tbFk8C3eYMw72NT5+WSu8/+/ZJBw86XxcVmfMrLjb/PnOmbH2Li31fBl+3aclch0eOSL/9Vn654mKzjuXVoaSk4ve9cR3v2DHvZU+e9LzPeZuvP20RLoqKAu9Xzif+7tORbPv2kM3ar+Bx+PBhFRcXKz4+3m14fHy8CgoKPI5TUFCghIQEn8sHzbhxUp06nt+LipJiYsy/L7lEio6WLrpIqlZN6tlTqlLFfC8mRoqLM/+Ojjbft8a3xMZKDRqY70dFSdWre55ndLRznq7jXnaZOTw62qxvw4aex61a1X3ernWpWlVKTzefS8/DUqWKc7msca1l9Mbbgb9WLeffcXFmvZs1M59dl79WLfeyFtf5W8vl6pJLpBo1pBtu8F630tMpPdxq77g4z8sRFWWuc0uVKp7Lx8WZy1arlrm+OnQoO70qVaS6dc3hrttWrVpmPcpTp4657lxVreo+j2rVzANJnTpmfR58UMrMlFz3y7p1zfnFxDjXs7WcVhtHR0sXXlh+faxyruuxomXwNL43MTFSzZrmIy7Oc9moqLL1tNqjenXnPlkea3+WzPVmbZ+lt1FPda9oea391VO9y9ufJOd25rqOXMf3xNM+4gvX6UVHm32bp/o1bGi2S6BB39/t42xU9oeRUH+4sVOLFiGb9fnzv1o6dpRSUswNKyVF6tdPyskxn9etk9asMf/+6CPzk9hnn0l33y09+6y0erX53po10ldfmX+vXWt+ojMMqW9f53y6dDH/J0FxsfnJqLDQnE/jxlJSklkmKkrq08f5KTQvz3y/qMj81N+woTn9X3+VEhPN8vHx5nNmpjntXr3MzqpvX3Mab7xhTmPuXHM6K1aYZSTzoG09Dx5slps92/wElJdnLmdenjkdq05WvVzbae5c8zkz05x3Zqb5+ssvneOcPGnWe+dO87mw0CwTFWWeMbjjDnN68fHmIyfHbEdrXrNmOeeZmWk+2rc3l+mHH5ztHB9vvmfVLS/PfTqlhxcXS/fcYy7fzTdL9eqZ04mNNbeHdeukQ4fMZWjc2CzXsKG5PHPnSvXrm+W++spctmPHpFOnpA0bzPVsreOcHLNdDxwwh5eUOLetL78013uVKmZZ1zZ84QVz2LRp5rrLzHSu76Iis35WezzzjNSokVn25Elz3S9fbj6sZT5wwKzjmTPm/2mIiTGnv26dc9vp00c6fNjZZtY24Lru8/LMtrPWozWep3FeeMH9QFavnnP8nBzzdWqqWc5qzzVrpDvvNNvy0kvNg5bVNklJ5jb7/PNmPV1PFFtnQu66y2zv+vWd+0hysrleJSkhwZzWLbc499dTp5zb56+/Sv/6l3ObSUkxp9WsmTn+tdea+2K/fuZwS40a5nJY/UJxcdk2sda9Nc3UVOeBLSHBXJerV5tlevd2Tj8uznxv3TpzetHR5rCUFHPaRUXubR0XZwY3yVnG2qZq1DDn2bixdN115t/165vL9P775j5RpYqzjikp5vBff3XOPzPTuf1Yy5aXZw5PTXU+W9vrTTc5/87Jca+DtV9a06xf3/zbdTuytn1r+ynvgoHrvpmXZz6sZbHWjzXPwYOd669+fff1Ze2/c+c6h7sut2t/Unof8bYsrtPPy3OuX6vvsuZbr565zWZmmnWuV88sExPj3AcGDzbft9az6zq3pKY65zl4sHt/7/pBt149s8y0aQoZf67fnDp1ymjVqpXx+eefuw0fOXKkMWTIEI/jpKenG3PmzHEbNnXqVKNHjx4eywf1Wy2+3N1rlbF+nriiO6Nd71gvfcOdpzuby/v2TGamea3e9Vsnnm4i9DT98u4qL+/O+PLu5PZWX1/r5G06/t497jq/in7voKLl8eVOb293zFfGD/F4azvXO9wrWtZA26/0evJ1/XmbTmnebj6taPzS9bGevf2wlrdvJZT3bQlfvx1i1cf1hmLXG4Ndb/j19G0Lb+1irVtP69R63/oqbUqKc7rl1Tsvz/2HyVxvAra+0ltR+7kO87Zd+dJvevrJ94puSPb122YVKV1vb+ulXz/z/pwaNXzbh13bt7x9zddl8WV6FfXZpedl3bvo+rPp3r61mZTk2ze8zkJQby4dO3as43VxcbHRuXPncm8uffjhh92G3XPPPfbfXOora8VZ32YpHTw8beSedtZAvpplGIF/9SkY4/XrV/7vm/jaFpVVV1/GD+av9AXzp4etZXL9+vHZHAh8Kedt/bkexHzh6SukvsyvIlZHHR1tPnv6arav0/SnLV3LVhRk/Tkw+Xow8XW63upf0fqoKAAGsl9WNE1P6+5s+wFv0/HWfoHM72z71dLL7e1DQHkfIP2pn7f6Bno88lNQv06blJRkLFmyxNi5c6cxevRoo2PHjkZ+fr5hGIbx1FNPGS+//LKj/IYNG4zWrVsbs2fPNnbu3GlMmzYtNF+n9VUgnURl7sSVtTNWBtdPfp4+MfnaFpWtvJ3obA/W5fE11JxNu3j61B5onSvqbLzV09/fXvE1kAXa8Xv6VO/vNP3pkCtaB4GuI386/7NZ9xXNx6aDkJuKAp4d87OTryG+dBsEq942tUdQf0Bs3rx5xo033mgkJiYaWVlZxqZNmxzvDRgwwBg1apRb+Y8//ti4+eabjcTEROO2224L7Q+IhbtQdAq+Cpe6+Xpq05WddQ/0TFBlsyv8VvZZJjsDbSDTDjRA2HUwrGg+oT4oh1s9gqGyzkpGGH4yPRQq41SlP9ei7RapO0npMzd2zC8UZ4LOFeHyibwyyoZLWAdsQPAIhbPpZM6XDioUn17DsW1DFUYiIQSFoo7B2kbCqb3D7cwRzjm+Hr89fHkcARs2zP3ZrnEjyZQp0qJF5t8LFtgz7XBs22C2QzjO1x8dO9pft2BtI6FYFm9Cse8BHhA8KtPZdDLh1EEFUzBDgLdph2PbhioMhWMICwfhuI1UtlDse4AHUYYRXr8Rm5GRIUnKzc0NcU0AAGFt/XrzbMuwYWZ4REj5evzmjAcAIDJxiScinT8/mY7gWr9e6t9fevtt83n9+lDX6NxhtS1tCrgbNsz8uXku8UQUznigclifPP79b/P/zUh8AqksfKoDPDsf7s05BxE8UDmsTxyZmdKnn577n0DsvLbMjXsAziEED1QO108eAweGti52sPMsBJ/qAJxDuMcj0nH9PzS4tgwAAeGMR6R77jnps8+kgwel5ctDXZvzB2chACAgnPEAAAC24YxHpBs3ToqP55Q/ACAicMYj0lmn/MPlV/u45wSILOyzsBlnPFC5+M0JILKwz8JmBA9ULn5zAogs7LOwGcEDlYtvewCRhX0WNuMeDwAAYBuCBwAAsA3Bw8Kd3QAABB33eFi4sxsAgKAjeFi4sxsAgKAjeFi4sxsAgKDjHg8AAGAbggcAALANwQMAANiG4AEAAGxD8AAAALYheAAAANsQPAAAgG0IHgAAwDYEDwAAYBuCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtgm7/0574MABFRcXKyMjI9RVAQAAPtq3b59iYmIqLBd2ZzyqVaumKlXCLg8BAIByVKlSRdWqVauwXJRhGIYN9QEAAAi/Mx4AAODcRfAAAAC2IXgAAADbEDwAAIBtCB4AAMA2ERs8fvnlF40YMULXXXed2rRpox49eui///2v1/Jr165Vy5Ytyzzy8/Pdyi1YsEBdu3ZVcnKyevXqpW+//TbYixLWgtHO06dPL/N+9+7d7VicsOVvO0tSUVGRpkyZoi5duigpKUldu3bVe++951bmk08+Uffu3ZWcnKwePXpo5cqVwVyMsBeMdl6yZEmZ7Tk5OTnYixL2/G3rp59+2mPfcdttt7mVo492F4x2DnYfHZE/mHHkyBH17dtX1113nd58801ddNFF+vHHH1WnTp0Kx12+fLlq1arleB0fH+/4++OPP1Z2drbGjBmjtm3bKicnR4MGDdLy5cvdyp0vgtXOktS8eXPNmTPH8dqXH505VwXazk8++aQOHjyo8ePH6/LLL1d+fr5KSkoc73/zzTcaPny4/vjHP6pLly768MMP9dhjj2nJkiVq0aJFsBcr7ASrnSWpVq1aWr58ueN1VFRUUJYhUgTS1s8++6yGDx/ueF1cXKw777zT7YBHH+0uWO0sBbmPNiLQpEmTjL59+/o1zpo1a4wWLVoYR44c8VomKyvLGDNmjON1cXGxkZaWZsycOTPgukayYLXztGnTjDvuuONsq3fOCKSdV65caVx99dXG4cOHvZZ58sknjYceeshtWK9evYzRo0cHUs2IF6x2Xrx4sXH11VefZe3OLYG0dWmff/650bJlS2PPnj2OYfTR7oLVzsHuoyPyUssXX3yhpKQkPfHEE0pNTVXPnj21aNEin8bt2bOn0tLSdP/992vDhg2O4UVFRdqyZYs6derkGBYdHa1OnTpp48aNlb4MkSAY7Wz58ccflZaWpoyMDA0fPlx79+6t7OpHjEDa2Rpn1qxZ6ty5szIzM/Xiiy/qt99+c5TZtGmTUlNT3cZLS0vTpk2bgrEYYS9Y7SxJhYWF6tKli9LT0/XII49ox44dwVyUsHc2fYflvffeU6dOndSgQQNJ9NGeBKOdLcHsoyMyeOzevVvvvvuumjRpotmzZ6tv374aN26c3n//fa/j1K1bV2PGjNG0adM0bdo01atXTwMHDtSWLVskSYcPH1ZxcXGZ03Xx8fEqKCgI6vKEq2C0syS1adNG2dnZmjVrll544QX9/PPP6t+/v44fP27HYoWdQNp59+7d2rBhg3bs2KEZM2boT3/6kz799FONGTPGUaagoEAJCQlu47E9V347X3HFFZowYYJee+01TZo0SYZhqE+fPtq/f78dixWWAmlrV7/88ov+9a9/KSsryzGMPrqsYLSzZEMfHbRzKUGUmJho3HPPPW7D/u///s/o3bu3X9Pp37+/MWLECMMwDGP//v1GixYtjG+++catzIsvvmhkZWWdXYUjVDDa2ZMjR44YHTp0MBYtWhRQPSNdIO18//33G8nJycbRo0cdwz799FOjZcuWxsmTJx3T/fDDD93Gmz9/vpGamlqJtY8cwWrn0oqKioxu3boZU6ZMqZR6R6Kz7Ttef/1149prrzVOnTrlGEYfXVYw2tmTyu6jI/KMR926ddWsWTO3YU2bNvX7VFBycrJ++uknSdJFF12kmJgYHTx40K3MwYMHy3xqPF8Eo509ueCCC9SkSZNyy5zLAmnnunXr6tJLL1Xt2rUdw5o1aybDMByftBMSEsp8EmR7rvx2Lq1q1apq1arVebs9S2fXdxiGocWLF+vOO+9UbGysYzh9dFnBaGdPKruPjsjg0aFDB33//fduw3744Ycy16gqsnXrVtWtW1eSFBsbq8TERK1evdrxfklJiVavXq327duffaUjUDDa2ZMTJ05o9+7d5ZY5lwXSzh06dNCBAwd04sQJx7Dvv/9e0dHRqlevniSpXbt2WrNmjdt4q1atUrt27Sqv8hEkWO1cWnFxsbZv337ebs/S2fUd69at048//ljm9D99dFnBaGdPKr2PrpTzJjb7z3/+Y7Ru3dr461//avzwww/GBx98YLRt29b4xz/+4Sjz8ssvG0899ZTj9Zw5c4zPP//c+OGHH4xt27YZ48aNM6666ipj1apVjjIfffSRkZSUZCxZssTYuXOnMXr0aKNjx45Gfn6+rcsXLoLVzhMnTjTWrl1r7N6929iwYYPx+9//3rjuuuuMgwcP2rp84SKQdj5+/Lhxww03GH/4wx+MHTt2GOvWrTNuvvlm49lnn3WU2bBhg9G6dWtj9uzZxs6dO41p06YZiYmJxrZt22xdvnARrHaePn268dVXXxk//fSTsXnzZmPYsGFGcnKysWPHDluXL5wE0taWESNGGL169fI4Xfpod8Fq52D30REZPAzDML744gvj9ttvN5KSkozu3bsbf/vb39zeHzVqlDFgwADH6zfeeMPo1q2bkZycbFx77bXGgAEDjNWrV5eZ7rx584wbb7zRSExMNLKysoxNmzYFfVnCWTDaeejQocb1119vJCYmGp07dzaGDh1q/Pjjj7YsT7jyt50NwzB27txp/P73vzfatGlj3HDDDUZ2dnaZ+w4+/vhj4+abbzYSExON2267zVixYkXQlyWcBaOdx48f7+gzOnXqZAwePNjYsmWLLcsTzgJp66NHjxpt2rQpU9YVfbS7YLRzsPvoKMMwjMo5dwIAAFC+iLzHAwAARCaCBwAAsA3BAwAA2IbgAQAAbEPwAAAAtiF4AAAA2xA8AACAbQgeAADANgQPAABgG4IHAACwDcEDAADY5v8DYbuq/GB2OTsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    feat = feat_mean + 3 * feat_std\n",
    "    feat_inv = feat[unq_inv]\n",
    "    mask = xyz[:, -1] > feat_inv\n",
    "\n",
    "rand_x, rand_y = np.random.randint(0, voxel_dim[0]), np.random.randint(0, voxel_dim[1])\n",
    "vox_mask = (voxel_index[:, 0] == rand_x) & (voxel_index[:, 1] == rand_y)\n",
    "\n",
    "if vox_mask.sum() == 0:\n",
    "    print(f\"No points in voxel ({rand_x}, {rand_y})\")\n",
    "else:\n",
    "    z_max = feat[unq_inv[vox_mask][0]]\n",
    "\n",
    "    plt.figure()\n",
    "    plt.style.use('seaborn-white')\n",
    "    # plt.axis('equal')\n",
    "    # plt.scatter(xyz[::500, 0], xyz[::500, 1], c='dimgrey', s=0.5)\n",
    "    plt.scatter(xyz[vox_mask, 0], xyz[vox_mask, 2], c='r', s=0.5)\n",
    "    plt.plot([xyz[vox_mask, 0].min(), xyz[vox_mask, 0].max()], [z_max, z_max], 'b--')\n",
    "    plt.title(f'Voxel Index: ({rand_x}, {rand_y})')\n",
    "    # plt.xlim(-9, 7.5)\n",
    "    # plt.ylim(-8, 8)\n",
    "    # plt.ylim(-7, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.3 Fly Around Certain Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "with torch.no_grad():\n",
    "    gaussians._scaling = org_scaling.clone()\n",
    "    gaussians._opacity = org_opacity.clone()\n",
    "    gaussians._scaling[mask] = -1e6\n",
    "    gaussians._opacity[mask] = -1e6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_70388/459775162.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    214\u001b[0m     \u001b[0;31m# for matrix city, z_dim=2, otherwise z_dim=1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    215\u001b[0m     \u001b[0mviewpoint_cam\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mloadCamV4\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mposes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxyz\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxyz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mangle\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mangle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 216\u001b[0;31m     \u001b[0mimg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrender\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mviewpoint_cam\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgaussians\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbackground\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"render\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    217\u001b[0m     \u001b[0mimg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mimg\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m255\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclamp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m255\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0muint8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpermute\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcpu\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    218\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/python_workspace/3DGS/gaussian_renderer/__init__.py\u001b[0m in \u001b[0;36mrender\u001b[0;34m(viewpoint_camera, pc, pipe, bg_color, scaling_modifier, override_color)\u001b[0m\n\u001b[1;32m     93\u001b[0m         \u001b[0mscales\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscales\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m         \u001b[0mrotations\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrotations\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 95\u001b[0;31m         cov3D_precomp = cov3D_precomp)\n\u001b[0m\u001b[1;32m     96\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     97\u001b[0m     \u001b[0;31m# Those Gaussians that were frustum culled or had a radius of 0 were not visible.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/gaussian_splatting/lib/python3.7/site-packages/torch/nn/modules/module.py\u001b[0m in \u001b[0;36m_call_impl\u001b[0;34m(self, *input, **kwargs)\u001b[0m\n\u001b[1;32m   1128\u001b[0m         if not (self._backward_hooks or self._forward_hooks or self._forward_pre_hooks or _global_backward_hooks\n\u001b[1;32m   1129\u001b[0m                 or _global_forward_hooks or _global_forward_pre_hooks):\n\u001b[0;32m-> 1130\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mforward_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1131\u001b[0m         \u001b[0;31m# Do not call functions when jit is used\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1132\u001b[0m         \u001b[0mfull_backward_hooks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnon_full_backward_hooks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/gaussian_splatting/lib/python3.7/site-packages/diff_gaussian_rasterization/__init__.py\u001b[0m in \u001b[0;36mforward\u001b[0;34m(self, means3D, means2D, opacities, shs, colors_precomp, scales, rotations, cov3D_precomp)\u001b[0m\n\u001b[1;32m    217\u001b[0m             \u001b[0mrotations\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    218\u001b[0m             \u001b[0mcov3D_precomp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 219\u001b[0;31m             \u001b[0mraster_settings\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    220\u001b[0m         )\n\u001b[1;32m    221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/gaussian_splatting/lib/python3.7/site-packages/diff_gaussian_rasterization/__init__.py\u001b[0m in \u001b[0;36mrasterize_gaussians\u001b[0;34m(means3D, means2D, sh, colors_precomp, opacities, scales, rotations, cov3Ds_precomp, raster_settings)\u001b[0m\n\u001b[1;32m     39\u001b[0m         \u001b[0mrotations\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m         \u001b[0mcov3Ds_precomp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m         \u001b[0mraster_settings\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m     )\n\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/gaussian_splatting/lib/python3.7/site-packages/diff_gaussian_rasterization/__init__.py\u001b[0m in \u001b[0;36mforward\u001b[0;34m(ctx, means3D, means2D, sh, colors_precomp, opacities, scales, rotations, cov3Ds_precomp, raster_settings)\u001b[0m\n\u001b[1;32m     90\u001b[0m                 \u001b[0;32mraise\u001b[0m \u001b[0mex\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     91\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 92\u001b[0;31m             \u001b[0mnum_rendered\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradii\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgeomBuffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbinningBuffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimgBuffer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_C\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrasterize_gaussians\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m         \u001b[0;31m# Keep relevant tensors for backward\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# find corner poses\n",
    "region_center1 = np.array([3.2, 2.8, 2.5])\n",
    "region_center2 = np.array([1.4, -2.75, 5.0])\n",
    "radius = 2\n",
    "\n",
    "poses = [views[i] for i in simplified_hull]\n",
    "\n",
    "video_path = os.path.join(lp.model_path, filename, \"ours_lod_video\")\n",
    "makedirs(video_path, exist_ok=True)\n",
    "\n",
    "# second parameter appoints number of points interpolated between each consecutive points to draw the curve\n",
    "xyz_pre = None\n",
    "\n",
    "fontsize = 20\n",
    "fig = plt.Figure(figsize=(12, 9), dpi=100)\n",
    "fig.set_tight_layout(True)\n",
    "canvas = FigureCanvasAgg(fig)\n",
    "ax = fig.gca()\n",
    "gs_xyz = gaussians.get_xyz.cpu().detach().numpy()\n",
    "ax.scatter(gs_xyz[::200, 0], gs_xyz[::200, dim], c='dimgrey', s=0.5)\n",
    "# ax.plot(cam_center_array[simplified_hull,0], cam_center_array[simplified_hull, dim], 'r--^',lw=2)\n",
    "ax.axis('equal')\n",
    "ax.set_xlabel('x/100m', fontsize=fontsize)\n",
    "ax.set_ylabel('y/100m', fontsize=fontsize)\n",
    "ax.tick_params(axis='both', which='major', labelsize=fontsize)\n",
    "ax.set_xlim(np.min(cam_center_array[:, 0]), np.max(cam_center_array[:, 0]))\n",
    "ax.set_ylim(np.min(cam_center_array[:, dim]), np.max(cam_center_array[:, dim]))\n",
    "\n",
    "# First Stage, fly around certain point\n",
    "\n",
    "# step = 0.01\n",
    "# plt.style.use('seaborn-white')\n",
    "# for t in np.arange(0, 1+step, step):\n",
    "#     xyz = region_center2.copy()\n",
    "#     xyz[0] += (radius + 1.0 * np.cos(2 * np.pi * t)) * np.cos(2 * np.pi * t + np.pi / 2)\n",
    "#     xyz[dim] += (radius + 1.0 * np.cos(2 * np.pi * t)) * np.sin(2 * np.pi * t + np.pi / 2)\n",
    "\n",
    "#     xyz[2] += 1.0 * np.cos(2 * np.pi * t)\n",
    "\n",
    "#     if t == 0 or t == 0.25 or t == 0.5 or t == 0.75 or t == 1:\n",
    "#         print(f\"Current Position: {xyz}\")\n",
    "\n",
    "#     delta = region_center2 - xyz\n",
    "    \n",
    "#     if not np.any(delta):\n",
    "#         continue\n",
    "    \n",
    "#     ax.scatter(xyz[0], xyz[dim], s=100, c='r', marker='x', label='poses')\n",
    "#     if t == 0:\n",
    "#         ax.legend(fontsize=fontsize)\n",
    "#     # canvas.draw()\n",
    "#     # buf = canvas.buffer_rgba()\n",
    "#     # bev_map = np.asarray(buf)[..., :3]\n",
    "\n",
    "#     angle = np.arctan2(-delta[0], delta[dim])\n",
    "\n",
    "#     # for matrix city, z_dim=2, otherwise z_dim=1\n",
    "#     viewpoint_cam = loadCamV3(lp, idx, poses[0], 1.0, xyz=xyz, z_dim=2, yaw=angle)\n",
    "#     img = render(viewpoint_cam, gaussians, pp, background)[\"render\"]\n",
    "#     img = (img * 255).clamp(0, 255).to(torch.uint8).permute(1, 2, 0).cpu().numpy()\n",
    "\n",
    "#     # img = np.concatenate([img, bev_map], axis=1)\n",
    "    \n",
    "#     frames.append(img)\n",
    "\n",
    "#     xyz_pre = xyz\n",
    "\n",
    "# Second Stage, fly around the whole scene\n",
    "\n",
    "pose_per_segment = 150\n",
    "pose_centers = np.array([\n",
    "    [3.2, 5.8, 3.5],\n",
    "    [1.2, 3.8, 2.5],\n",
    "    [3.2, 1.8, 1.8],\n",
    "    [5.2, 3.8, 2.5],\n",
    "    [3.2, 5.8, 3.5],\n",
    "    [-3.0, 4.1, 5.0],\n",
    "    [-7.37, -0.46, 5.0],\n",
    "    [-3.16, -4.7,  6.0],\n",
    "    [3.5, -5.74, 5.0],\n",
    "    [7, -2, 5.0],\n",
    "    [5.0, 1.0, 3.0],\n",
    "    [2.0, 0.5, 2.0],\n",
    "    [-3.0, 2.9, 2.0],\n",
    "    [-5.5, 1.2, 2.0],\n",
    "    [-4.0, -1, 3.5],\n",
    "    [1.0, 0.0, 6],\n",
    "    [3.5, -1.8, 5.5],\n",
    "    [1.0, -3.6, 4.5],\n",
    "    [0.5, 0, 6],\n",
    "    [0.5, 0, 10],\n",
    "    [0.5, 0, 15],\n",
    "])\n",
    "\n",
    "interp_xyz, acc_pts = evaluate_bezier_v3(pose_centers, pose_per_segment)\n",
    "\n",
    "Rt = np.zeros((4, 4))\n",
    "Rt[:3, :3] = poses[0].R.transpose()\n",
    "Rt[:3, 3] = poses[0].T\n",
    "Rt[3, 3] = 1.0\n",
    "\n",
    "C2W = np.linalg.inv(Rt)\n",
    "ref_angle = np.array(mat2euler(C2W[:3, :3]))\n",
    "\n",
    "key_point1, interval1 = acc_pts[3], 10\n",
    "key_point2, interval2 = acc_pts[9], 50\n",
    "key_point3, interval3 = acc_pts[14], 50\n",
    "key_point4, interval4 = acc_pts[17], 50\n",
    "# key_point2, interval2 = acc_pts[2], 50\n",
    "# key_point3, interval3 = acc_pts[7], 50\n",
    "# key_point4, interval4 = acc_pts[10], 50\n",
    "xyz_pre = None\n",
    "\n",
    "\n",
    "xyz_list, angle_list = [], []\n",
    "\n",
    "for t in range(len(interp_xyz)):\n",
    "    xyz = interp_xyz[t]\n",
    "    angle = ref_angle.copy()\n",
    "\n",
    "    if xyz_pre is not None and not np.any(xyz - xyz_pre):\n",
    "        continue\n",
    "\n",
    "    if t < key_point1 - 10:\n",
    "        # stage 1\n",
    "        delta = region_center1 - xyz\n",
    "        angle[2] = np.arctan2(delta[0], -delta[dim]) + np.pi\n",
    "        \n",
    "    elif t < key_point1:\n",
    "        # transition\n",
    "        delta = region_center1 - xyz\n",
    "        angle2 = np.arctan2(delta[0], -delta[dim]) + np.pi\n",
    "\n",
    "        delta = xyz - xyz_pre\n",
    "        angle1 = np.arctan2(delta[dim], delta[0])\n",
    "\n",
    "        angle[2] = (1 - (t - key_point1 + interval1 + 1) / interval1) * angle2 + (t - key_point1 + interval1 + 1) / interval1 * angle1\n",
    "\n",
    "        # print(f\"circle: {angle2}, org: {angle1}, current: {angle}\")\n",
    "    \n",
    "    elif t < key_point2 - interval2:\n",
    "        # stage 2\n",
    "\n",
    "        if xyz_pre is None:\n",
    "            delta = interp_xyz[t + 1] - interp_xyz[t]\n",
    "        else:\n",
    "            delta = xyz - xyz_pre\n",
    "        \n",
    "        # delta = xyz - xyz_pre\n",
    "        angle[2] = np.arctan2(delta[dim], delta[0])\n",
    "\n",
    "    elif t < key_point2:\n",
    "        delta = xyz - xyz_pre\n",
    "        angle2 = np.arctan2(delta[dim], delta[0])\n",
    "        angle1 = np.arctan2(-delta[0], delta[dim])\n",
    "        \n",
    "        angle[2] = (1 - (t - key_point2 + interval2 + 1) / interval2) * angle2 + (t - key_point2 + interval2 + 1) / interval2 * angle1\n",
    "\n",
    "    elif t < key_point3 - interval3:\n",
    "        # stage 2\n",
    "        delta = xyz - xyz_pre\n",
    "        angle[2] = np.arctan2(-delta[0], delta[dim])\n",
    "\n",
    "    elif t < key_point3:\n",
    "        # transition\n",
    "        \n",
    "        delta = xyz - xyz_pre\n",
    "        angle2 = np.arctan2(-delta[0], delta[dim])\n",
    "\n",
    "        delta = region_center2 - xyz\n",
    "        angle1 = np.arctan2(delta[0], -delta[dim]) - np.pi\n",
    "\n",
    "        angle[2] = (1 - (t - key_point3 + interval3 + 1) / interval3) * angle2 + (t - key_point3 + interval3 + 1) / interval3 * angle1\n",
    "    \n",
    "    elif t < key_point4 - interval4:\n",
    "        # stage 3\n",
    "        delta = region_center2 - xyz\n",
    "        angle[2] = np.arctan2(delta[0], -delta[dim]) - np.pi\n",
    "    \n",
    "    elif t < key_point4:\n",
    "        # transition\n",
    "        delta = region_center2 - xyz\n",
    "        angle[2] = np.arctan2(delta[0], -delta[dim]) - np.pi\n",
    "\n",
    "        angle[0] = (1 - (t - key_point4 + interval4 + 1) / interval4) * angle[0] + (t - key_point4 + interval4 + 1) / interval4 * (-np.pi)\n",
    "    \n",
    "    else:\n",
    "        # stage 3\n",
    "        delta = region_center2 - xyz\n",
    "        angle[2] = np.arctan2(delta[0], -delta[dim]) - np.pi\n",
    "        angle[0] = -np.pi\n",
    "\n",
    "    xyz_list.append(xyz)\n",
    "    angle_list.append(angle)\n",
    "\n",
    "    xyz_pre = xyz\n",
    "\n",
    "angle_list = np_move_avg_v2(np.unwrap(np.array(angle_list), axis=0), 10, mode='same').tolist()\n",
    "\n",
    "frames = []\n",
    "for t in range(len(xyz_list)):\n",
    "    xyz = xyz_list[t]\n",
    "    angle = angle_list[t]\n",
    "\n",
    "    ax.scatter(xyz[0], xyz[dim], s=100, c='r', marker='x', label='poses')\n",
    "\n",
    "    if t == 0:\n",
    "        ax.legend(fontsize=fontsize)\n",
    "    \n",
    "    # canvas.draw()\n",
    "    # buf = canvas.buffer_rgba()\n",
    "    # bev_map = np.asarray(buf)[..., :3]\n",
    "\n",
    "    # for matrix city, z_dim=2, otherwise z_dim=1\n",
    "    viewpoint_cam = loadCamV4(lp, idx, poses[0], 1.0, xyz=xyz, angle=angle)\n",
    "    img = render(viewpoint_cam, gaussians, pp, background)[\"render\"]\n",
    "    img = (img * 255).clamp(0, 255).to(torch.uint8).permute(1, 2, 0).cpu().numpy()\n",
    "\n",
    "    # img = np.concatenate([img, bev_map], axis=1)\n",
    "    \n",
    "    frames.append(img) \n",
    "\n",
    "# canvas.draw()\n",
    "# buf = canvas.buffer_rgba()\n",
    "# bev_map = np.asarray(buf)[..., :3]\n",
    "# plt.imshow(bev_map)\n",
    "# plt.axis('off')\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(np.array(angle_list)[:, 2] * 180 / np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.5 Save Video"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "IMAGEIO FFMPEG_WRITER WARNING: input image is not divisible by macro_block_size=16, resizing from (1600, 900) to (1600, 912) to ensure video compatibility with most codecs and players. To prevent resizing, make your input image divisible by the macro_block_size or set the macro_block_size to 1 (risking incompatibility).\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Video saved to ../output/block_mc_aerial_block_all_lr_c36_loss_5_50_lr64_vq/block_all_test/ours_lod_video\n"
     ]
    }
   ],
   "source": [
    "video = imageio.get_writer(os.path.join(video_path, \"video.mp4\"), mode=\"I\", fps=30, codec=\"libx264\", bitrate=\"16M\", quality=10)\n",
    "for frame in frames:\n",
    "    video.append_data(frame)\n",
    "video.close()\n",
    "print(f'Video saved to {video_path}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gaussian_splatting",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
